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ABSTRACT 


An Interactive computer proRrnm (In FORTRAN) for the 
varin/iice-covarlnuce analyslH of VUll experiments Is pre»;:nted for use 
In experiment piannlnR, simulation studies and optimal design problems. 
The interactive mode la especially suited to these types of analyses 
providing ease of operation us well as savings in time and cost. The 
geodetic parameters Include baseline vector parameters and variations 
in polar motion and earth rotation. 

A discussion of the theory on which the program is based pro- 
vides an overview of the VLBI process emphasizing the areas of interest 
to geodesy. Special emphasis is placed on the problem of determining 
correlations between simultaneous observations from a network of 
stations. A model suitable for covariance analyses is presented. Sug- 
gestions towards developing optimal obsoi*vatlon schedules are Included. 


iv 




ACKNOWLEDGEMENTS 

I wish to express my gratitude to my advisor, Dr. Ivan I. 

Mueller for his guidance, encouragement, useful suggestions and 
patience throughout the course of my studies and the preparation of this 
report. I appreciate his efforts to provide me with the opportunity to 
meet with people working in the VLBI field. I am grateful to 
Dr. Richard H. Riipp for reading the report and offering valuable sug- 
gestions. 

Special thanks are due. to Dr. Douglas S, Robertson of the 
National Geodetic Survey and Dr. Chopo Ma of Goddard Space Flight Center 
for answering my numy questions pertaining to VLBI. I wish to thank 
Dr. Irwin I. Shapiro of M.I.T. for taking interest in my work and 
offering enlightening comments. I am grateful for the help provided me 
by Dr. Chopo Ma, Dr. Tom A. Clark, Mr. Jim Ryan, Mr. Bruce 5 . Schupler 
and Dr. Nancy Vandenbe.rg of GSFC during ray visit with them. I appreci- 
ate the use of the $77JUN26 and $78MAY17XC VLBI data sets. 

1 am very grateful for the financial support received as a 
Graduate Research Associate at The Ohio State University. I am also 
indebted to the Instruction and Research Center at The Ohio State 
University for providing extensive computer support. 

Thanks are due to Mr, KitIcos Pavlis and Mr. Lenny Krleg for 
reading through the original manuscript and making useful comments, to 


V 




Ms. Irene Tesfal and Ms. Sandy Smith for preliminary typing and to 
Ms. Linda Wright for the typing of the final manuscript. 

I wish to thank all my friends at the Department of Geodetic 
Science for making my stay here till now a fruitful and pleasant one. 


TABMl OF CONTENTS 


Page 

DEDICATION il 

PREFACE til 

ABSTRACT iv 

ACKNOWLEDGEMENTS v 

LIST OF TABLES ix 

LIST OF FIGURES x 

1. INTRODUCTION 1 

1.1 Background . . 1 

1.2 Purpose of the Report 5 

1.3 Organization and Scope 7 

2. THE .VffiASUREMENT PROCESS 9 

2.1 Introduction 9 

2.2 Observables 9 

2.2.1 Basic Observables 9 

2.2.2 Geometric Observables 11 

2.2.3 Measured Observables . 15 

2.3 Data Acquisition and Observable Estimation 17 

2.4 Deviations from the Geometric Model ........ 22 

3. MATHEMATICAL MODELS 26 

3.1 Introduction 26 

3.2 Least Squares Adjustment Mathematical Models .... 27 

3.2.1 Introduction ....... 27 

3.2.2 Coordinate Systems Definition 27 

3.2.3 Time Delay Model 29 

3.2.4 Time Delay Rate Model 40 

3.2.5 Adjustment Algorithm 43 

3.2.6 Weighting of Observations 47 

3.2.7 Model Refinements 51 


vii 


3.3 Singularity Problems ..... 53 

3.3.1 Introduction ....... 53 

3.3.2 Observation Singularities 53 

3.3.3 Critical Configurations 56 

3.4 Source Visibility Equation 59 

4. CONCLUSIONS AND FUTURE RESEARCH 61 

4.1 Summary and Conclusions 61 

4.2 Optimal Design Problems 63 

4.3 Observation Correlations 70 

APPENDIX A; VlJl COVARIANCE ANALYSIS INTERACTIVE PROGRAM (VIP). 75 

A.l Introduction 75 

A. 2 Job Control Language (JCL) 76 

A. 3 Explanatory Information 60 

A. 4 VIP Documented Listing 69 

APPENDIX B; VIP SAMPLE RUN 138 

B. l Introduction 138 

B.2 A Typical Interactive Session • 141 

B.3 Post-Session Output 173 

REFERENCES 165 


vili 


LIST OF TABLES 


Tabla Paga 

A.l Comparing Simultaneous 3-Station Observations Using 
All 3 Baselines with 2-Baseline Combinations 
(Diagonal Weight Matrix - Unweighted Mode) 72 

4.2 Con^)arlng Simultaneous 3-Statlon Observations Using 
All 3 Baselines with 2-Basellne Combinations 
(Diagonal Weight Matrix - Weighted Mode) 73 

A.l VIP File Allocation . 78 

A. 2 VIP Subroutine Index . , , 82 

A. 3 VIP Input Parameters . 85 

A. 4 VIP Program Options 87 


ix 


I 


LIST OF FIGURES 


Figure 

2.1 Data Acquisition and Processing Flow for VLSI 

2.2 VLSI Geometry on a Rotating (Earth) Platform 


2.3 Geometry of a Time Delay Observation 


2.4 Mark III Field System 


A.l VIP CLIST 


A. 2 Program Flow 




1. IKrRODUGTXOK 


"A marriage of convenience has been consummated between 
the disparate fields of geodesy and radio astronon^. The 
radio technique of very-long-basellne interferometry 
(VLBI) promises to have a profom effect on studies of 
the Earth, Whether such promises will be fulfilled 
remains to be seen." 

[Counselman and Shapiro, 1978b] 

1.1 Background 

The application of VLBI to geodesy, geodynamics, and geo- 
physics Is an outgrowth of developments in the fields of radio astron- 
omy. In order to , >taln fine angular resolution In the study of the 
structure and size of extragalactlc radio sources, discovered by Jansky 
in the early 1930 's [Kraus, 1966], radio astronomers turned to inter- 
ferometry. In conventional interferometry, two (or more) radio anten- 
nas, acting as one impractlcally large single antenna (Increasingly 
finer angular resolution is roughly proportional to antenna diameter) 
are connected by cables whereby signals received from a radio source are 
compared instantaneously. The maximum separation in this mode Is 10-20 
km. With the development of very stable frequency standards and 
wide-band tape recorders, the real-time link between the antennas could 
be eliminated. Thus, the concept of very-long-basellne interferometry 
where antenna separation of thousands of kilometers is possible. In 
this mode, the received radio signals are tape recorded at each site 


1 


and cro88~corralatad Xatar at a cantral procasaing facility to recover 
the VLSI "obaervablea" deccribed in Chapter 2. 

VLB! meaauretnents , besides their radio-astronomy applications, 
supply information about baseline components and distances, radio 
source coordinates, polar motion, UTl, precession and nutation, and 
solid earth tides. However, except fov baseline distances and 
radio-source declinations, all the remaining geodetically relevant para- 
meters are non-estlmable unless defined as variations (in polar motion, 
etc.) as will be explained in Chapter 3. Anticipated observational accu- 
racies, resulting from increasingly better instrumentation, and improved 
mathematical models should make the monicoring of a global tectonic 
plate motions, continental drift and crustal deformations feasible. 

Tlie astronomic applications of VLBl include the possibility of 
classifying a number of radio sources as fixed, in order to define an 
absolute extra-galactic coordinate system, A list of such sources has 
been proposed in [Elsroore and Ryle, 1976J. There are presently several 
hundred known radio sources- With the Mark I recording system less 
than 20 sources were of sufficient strength for geodetic applications. 
With the state-of-the-art Mark III recording and processing system (see 
[Ma, 1978; Clark, 1979a] for .. «»*acriptlon) , this number should be 
increased considerably, thereby improving the distribution of sources 
in the sky. 

The first VLBI experiments were conducted in the late 1960*8 
[Broten et al , , 1967; Bare et al,, 1967; Hinteregger et al., 1972]. 
Since then, many experiments have been performed by groups in the 
United States (the "East coast" and "West coast"), Canada and Europe, 
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too nunuirous for all to ba daacribad hart. At praaont, baaallna 
langtha of up to aavaral thousand kllomatare hava baan maasurad with a 
rapaatabillty of undar 5 cm [Shapiro, 1978] , A 1.24 km baaallna vector 
has baan datarmlnad with approxlmataly 5 urn rapaatabillty from VLSI 
obsarvatlons [Rogars at al. , 1978] . Maasuramants of tha sama basallna 
by convantlonal gaodetlc tachniquas comparad encouragingly well at the 
few millimeter level [Carter, 1979] . A 42~km baseline measured 
with portable VLBI antennas and conventional methods compared to within 
a decimeter In length [Nlell, 1979]. Variations In polar motion and 
UTl have been estimated at the decimeter and millisecond levels respec- 
tively . In a series of experiments, polar motion results agreed well 
with IPMS and Doppler results, but a systematic difference was detected 
with BIH values [Robertson et al., 1978]. Prom the same experiments, 
VLBI and BIH UTl results were found to differ by « .ut 0?002 rms. 
Fanselow et al. [1979] reported measurements of earth rotation para- 
meters at the OVOl accuracy level that compared well with lunar laser 
ranging (LURE) results. Although the accuracy of VT.BI parameter esti- 
mation is not yet completely clear, the above results are good indica- 
tions that "promises will be fulfilled." Recent VLBI-satellite laser 
Intercomparison experiments using the Mark III system will provide inde- 
pendent checks on accuracy and hopefully point to unmodelled systematic 
errors .■ At present, the primary limiting factors on accuracy are clock 
behavior and the propagation medium, particularly the wet component of 
the troposphere. Other factors include uncalibrated instrumental 
errors, : ladequate modelling of geophysical and relativistic effects, 
source structure and gravitational flexure of the larger telescopes. 
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Besides Its excellent angular resolution and impressive accu- 
racy, VLBI provides several other advantages. VLSI measurements are 
Independent of the Earth’s gravity field. In addition, since the 
anten, ^s receive microwave radiation, VLBI has practically all weather 
capabilities. On the other hand, the necessary equipment Is expensive 
and the availability of permanent antennas Is limited. The latter 
problem can be remedied by the use of portable antennas such as those of 
the Astronomical Interferometric Earth Survey (ARIES) system [MacDoran 
et al., 1978], An interesting list of radio Interferometry ’'advantages 
and disadvantages" as well as for Doppler, satellite laser ranging and 
lunar laser ranging techniques Is given in [Mulholland, 1978]. 

The next decade will see VLBI move into tho operational stage as 
the following examples illustrate. The NASA Geodynamics Program will 
concentrate on the detection of crustal movements by VLBI and satellite 
laser techniques. A Crustal Dynamics Project has been established at 
Goddard Space Flight Center for this purpose [NASA, 1979a]. The Polar 
Motion Analysis by Radio Interferometric Surveying (Polaris) project is 
planned for the early 1980's [Carter, 1978]. Its goal is to establish 
and operate a three-scat Ion VLBI network to monitor polar motion and 
earth rotation on a regular basis. It la anticipated that small port- 
able interferometer terminals, receiving signals from Global Positioning 
System (GPS) satellites, will yield several millimeter accuracy for 
baseline lengths up to several hundred kilometers. These systems, oper- 
ating on the same basic principles of VLBI as described in this thesis, 
are now being developed. They include Miniature Interferometer Terminals 
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for Earth Surveying (MITES) [Counaelman and Shapiro, 19'.> :a] and Satellite 
Emieaion Radio Interferometric Earth Surveying (SERIES) [MacDoran, 1978]. 

For a detailed history of VLB I development as well as an exten- 
sive bibliography, see [Benjauthrit, 1978a and b]. A list of the various 
agencies participating In VLBI development and a description of several 
of their experiments are found in [Campbell, 1979] . Fanselow [1978] 
gives a summary of completed as well as current VLBI programs. The 
proceedings of the Radio Interferometry Techniques for Geodesy con- 
ference contains the most up-to-date description of the present status 
of VLBI [NASA, In presn]. Other valuable references, especially for 
geodetic applications, include [Thomas, 1972a and b, 1973; Whitney, 

1974; Robertson, 1975; Dermanls, 1977; Ma, 1978; Shapiro, 1978]. 

1.2 Purpose of the Report 

A VLBI covariance analysis Interactive Program (VIP) is pti,.'- 
sented for use in simulating and planning VLBI experiments. An explana- 
tion of the theory and mathematical models on which this program is 
based is intended to provide an overview of VLBI for those Interested in 
applying the VLBI technique to geodetic activities. 

VIP provides an upper limit on accuracy attainable for the VIP 
parameter set given the planned station configuration and source 
schedule of a particular experiment and the a priori noise estimates of 
delay and delay rate measurements. Only random errors are assumed and 
there is no provision for systematic effects except for a simple 
two-term polynomial to model errant clock behavior. Therefore, it is 
not expected that this type of analysis will reflect the actual 
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performance of a particular experiment which may be several times worse 
than the a priori numbers Indicate. Nevertheless, a covariance analy- 
sis is useful in comparing the relative effects of different station 
locations and observation schedules. Ultimately, the geometrical 
strength of a given experiment is of primary importance in optimal para- 
meter estimation. 

llie choice of parameter set was influenced by studies of dif- 
ferent observation schedules for the Polaris network mentioned above. 
Therefore, the main emphasis is on estimation of earth orientation para- 
meters including variations in polar motion and earth rotation as well 
us on baseline vector parameters. 

At its early stages of development, VIP was run in the batch 
mode. It was decided to modify the various routines to run in the 
interactive mode using the Time Sharing Option (TSO) and Tektronix ter- 
minals at the OSU Instruction and Research Computer Center (TRCC). In 
this mode, the user is able to simulate an experiment, view the results 
in real-time, and rerun through the program with the option of changing 
any or all of the initial input parameters. This process may be 
repeated as many times as desired with one loading of the program. 

Thus, the interactive mode is found to be ideal for this type of analy- 
sis, offering case and flexibility of operation as well as savings in 
time and cost. 

In all of the modern geodetic "space" systems, the geodesist has 
moved further away from the actual measurement process. In VLBI we are 
presented with a list of "observables," themselves estimated by a com- 
plex procedure requiring sophisticated instrumentation developed by 
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electrical engineers and radio astronomers. It Is Important to obtain 
familiarity with this measurement process (summarized In Chapter 2). 

With this background, the geodesist can address such problems as opti- 
mal experiment simulation and planning, development of Improved mathe- 
matical models, sound statistical analysis of data and correct 
adjustment philosophy, and, flna^j', can apply VLSI data to geodesy and 
Its related fields. Thesa probl m ' eas will be discussed and topics 
for future research presented. 

1.3 Organization and Scope 

Chapter 2 coveru the basic geometry of VLBI observations, the 
necessary Instrumentation, and explains the process by which the raw 
observed data Is transformed Into the "observables" of the least 
squares adjustment from which the geodetic parameters are estimated. 

In Chapter 3 the mathematical models used In VIP are described as well 
as possible model refinements. A sununary of VLBI estimable parameters 
Is Included. A model, suitable for covariance analyses, is presented 
for determining the correlations between simultaneous VLBI observations 
at a given epoch. Singularity problems arising from coordinate system 
definition, observability conditions and critical configurations are 
enumerated. An approach to observation schedule optimization is 
described in Chapter A. This last chapter also discusses the prob- 
lems related to obtaining correlations between simultaneous VLBI 
observations. Appendix A includes a documented listing of VIP plus 
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explanatory tables and figures. Appendix B contains the hard copy of a 
sample run as viewed on the Interactive screen. 
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2. THE MEASUREMENT PROCESS 


2.1 Introduction 

In this chapter the basic VLBI observables are described. 

First* their purely geometric Interpretations are presented followed by 
a discussion of the quantities that are actually measured. A brief 
description of the VLBI hardware is given* as well as the process by 
which the observables are estimated. This will be of a general nature 
only and the technical details may be found In the references. Expres- 
sions for the precision of the observables are included. Finally, sys- 
tematic errors that affect the measurement process and, thus* the 
estimation of geodetic parameters are summarized. Figure 2.1 (from 
[Fanselow* 1978]) Illustrates the measurement process and* therefore, 
the contents of this chapter. The parameter solution will be described 
in the next chapter. 

2.2 Observables 

2.2.1 Basic Observables 

A VLBI baseline consists of one antenna at each end* simul- 
taneously observing the random radio signals emitted from a compact 
extra-galactic source (e.g., a quasar). A particular segment of a wave- 
front will arrive at one antenna later than at the other as a result of 
the difference in path length to each station. Tills time delay is 
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Figure 2.1 





primarily a function of the location of the source In the extra-galactic 
frame and the baseline vector fixed to the rotating deformable earth. 

In Figure 2.2, f*om (Ma, 1978], we view the equatorial projection of a 
VLBI baseline at two different epochs. The path length difference Is 
given by cTi and ct^, respectively, Ti and Tjs are the time delays at the 
two epochs and c, the speed of light. As can be seen In the flguru, the 
time delay changes with time. Its rate of change Is called the time 
delay rate. 

Hie time delay and time delay rate contain the geodetlcally 
relevant Information. Any phenomena that affects these quantities can 
be theoretically parameterized In the mathematical model. The orienta- 
tion of the baseline with respect to the "inertial" frame is affected by 
polar motion and UTl variations. Therefore, the observables are sensi- 
tive to these changes although not to the absolute orientation of the 
baseline, as will be explained in the next chapter. The baseline vector 
is affected by solid earth tides and geodynamlc phenomena such as crus- 
tal motion. The source unit vector is affected by precession and nuta- 
tion. The estimable parameters will be defined in the next chapter, but 
it suffices to mention here that the observables are sensitive to these 
and other phenomena as well as to baseline vector and source coordir Ate 
parameters . 


2.2.2 Geometric Observables 

In this section, the geometric definition of the observables are 
presented under the assumption of perfect Instrumentation and of radio 
waves propagating in vacuum from a point source. The actual physical 
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conditions ar«, of course, quite different making the measured observ- 
ables vary considerably from their geometric counterparts. Some of 
these effects may be modelled better than others but all serve to com- 
plicate geodetic parameter estimation. They are described in section 
2.4. The measured observables, as will be seen in the next two sections 
are estimated by cross correlation of the tape recordings of the 
received signals. 

The basic geometry for a typical baseline Is shown in the 
figure below. 

^ SOURCE 



Figure 2.3. Geometry of a Time Delay Observation 

A certain segment of a wavefront arrives at antenna I at time ti 
and antenna 2 at time tj. The station 1 vector at time ti Is Xi(ti), 
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th 0 station 2 vector at tine tt» XgCts). The station position vectors 
are assumed for this discussion to be given in a geocentric Cartesian 
reference frame fixed with respect to the radio sources (assumed to be 
an inertial frame). From Figure 2.3 

Di 2 - -[XzCta) - X,(tx)] • 8 (2.2-1) 

where s is the unit vector in the direction of the source. The geo- 
metric time delay ls» therefore 

T " t2 “ t| =■ (2.2-2) 

g c 

In the Interval of time, T , station 2 has rotated by a small amount due 

s 

to the earth’s rotation. Since T is small (its maximum value is 

o 

approximately 0.02 see) we can write as a linear approximation 

X2(t2) - X2(tj) (2.2-3) 


From (2.2-1) - (2.2-3) 


T 

g 


c 


T_V2(tl) 

Bi I . 

c 
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where 

Bi2(ti) - X2(tj) - Xi(tj) 


(2.2-4) 


(2.2-5) 


is the instantaneous baseline vector at epoch tj and 

V2(ti) - - n X X2(ti) (2.2-6) 

the earth rotation vector (at ti). Notice that V 2 (tx)/c multiplied 
by the frequency of the received signal is the Doppler frequency shift. 
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In tht rnmalndar of the th«8ia» th« apaad of light will b« »*t 
to unity so that the time delay will be expressed in units of distance. 
From (2.2~4) the time delay is seen to be composed of two parts. 1110 
first term is the projection of the instantaneous baseline vector (at 
t|) in the direction of the source. The second term is the motion of 
station 2 during the wave transit. It is of small magnitude and can be 
accurately calculated based on a priori Information. Therefore, it can 
be neglected in developing the mathematical models in the next chapter. 
The time delay which is now In distance units will be expressed there as 


■“uk - • *k 

til tH 

where the subscript 1 refers to the 1 baseline, k to the k source 
th 

and j to the j epoch of observatioj . 

The time delay rate is then 


dB (t ) _ 

‘*ljk " ” dtj ' \ 

assuming that Sj^ ■* 0. 


( 2 . 2 - 8 ) 


2.2.3 Measured Observables 

The velocity of electromagnetic radiation passing through the 
atmosphere (a dispersive medium) can be divided into two categories, 
the group velocity and the phase velocity. Therefore, measurement of 
the difference in times of arrival may be of two types; the phase delay 
difference (called simply the phase delay) or the group delay difference 
(group delay) (Shapiro, 1978]. Theoretically, the phase delay could be 
calculated by dividing the phase difference of the recorded data streams 
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(call«d the fringt ph«s«) at a particular epoch by the (angular) fre- 
quency of the Incoming signal. However, the fringe phase Is asdilguous 
to some Integer multiple of 2ir, thereby Inflicting closely spaced ambi- 
guities on the phase delays which are difficult to resolve. The group 
delay, the derivative of the fringe phase with respect to angular fre- 
quency can be, theoretically, estimated unambiguously by measuring 
fringe phase over a wide band of frequencies. A simple example, based 
on a discussion by iMollnder, 1978], will Illustrate thesf points. 
Suppose that (|)^^ and are the fringe phases at frequencies f} and 
fj. Then 


(|>f^(t) - 2irfiT + 2iTm 
(j)^^(t) « 2TTf2T + 2irn 


(2.2-9) 


where T Is the time delay, 2im and 2Trn the ambiguities, m and n Inte- 
gers. If the uncertainty In the slope of fringe phase versus frequency 
is less than 2ir/(f2 -f j) then the ambiguities may be rescived and the 
time delay is given by 


T 


•I'e, - »f, 

2ir(fj-fi) 


( 2 . 2 - 10 ) 


Thus, fi and f2 must be spaced close enough so that the ambiguities may 
be resolved based on a priori information. A third frequency fs can 
then be spaced at an interval larger than f2 - f| because of the more 
accurate slope available from the previous determination. This proce- 
dure can be extended over several frequency bands, thus, the bandwidth 
synthesis technique [Rogers, 1970; Whitney et al,, 1976]. The wider the 
bandwidth, the more accurate the measurement of group delay. In the 
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Mark III isyatem, 28 narrow fraquancy banda, each 2 MHZ wide, era dia- 
trlbutad over a total of up to 400 MHZ [Shapiro, 1978]. 

Thua, the group delay la the maaaurad time delay. In practice, 
the do have andiigultlaa but thaae can be eliminated by 

examination of their reaiduala from an initial leaat aquarea adjuatment 
[Robertaon, 1975]. 

The fringe rate la the aecond, and leaa important, eatimated 
claervable. It ia the time derivative of the fringe phase. We will 
deal with the phase delay rate which la the fringe rate divided by the 
angular frequency. The phase delay rate la the measured time delay 
rate. One advantage of the phase delay rate (or the fringe rate) is 
that it can be determined unambiguously without resortl g to bandwidth 
synthesis, and therefore requires relatively simple equipment. How- 
ever, it suffers from several geometric disadvantages described in 
section 3.2.4 and is much less precise compared to the group delay. 

From this point on, we shall use Che terms delay and delay rate 
for the measured observables. 

2.3 Data Acquisition and Observable Estimation 

The Mark III field system (s<'u Figure 2.4 taken from [Ma, 1978]) 
is the state of the art in VLBI data acquisition hardware. This system, 
in conjunction with a radio antenna and environmental sensors, consists 
of basically a receiver, a frequency standard, a recorder and a phase 
calibrator. The entire system is run by the VLBI controller, an HP 1000 
mini-computer. Using schedule Input, the controller sets the receiver 
and recorder configurations, directs the telescope to a particular 
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Figure 2.4 Mark III Field System 























80 urc«i starts and stops the data drives, monitors the system's func- 
tions and logs all necessary Information [Ma, 1978]. 

Local oscillator signals, derived from the frequency standard, 
are mixed with the received radio-frequency signals. Several or all of 
the 28 possible channels are selected at all stations and are sampled, 
one channel per record of tape. The resulting intermediate-frequency 
signals are converted to video signals which are recorded on 
magnetic tape. For each tape record, the time epoch, derived from 
the frequency standard, is recorded. See [Whitney et al., 1976] for 
a detailed description of the system components. 

The phase calibration system Is used to reduce the dispersive 
effects of the instrumentation and to measure the timing cable length 
[Whitney et al., 1976; Thomas, 1978; Rogers, 1979]. 

At each station environmental sensors record temperature, 
humidity and pressure. A water vapor radiometer, if available, measures 
water vapor path delay [Claflln et al., 1978; Reach et al., 1979; Moran, 
1979]. 

The tapes from the participating stations are sent to a central 
processing facility for cross correlation. This Involves reconstructing 

A 

the radio-signal transmission process. A model delay T Is computed to a 
good approximation based on a priori Information. The data streams from 

A 

two tapes are offset by X and the signals are multiplied together. The 
theoretical cross correlation function is given by 

00 ^ y 

f X,(t)Xa[t + (x+'rt)]dt (2.3-1) 

^00 
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where X} and Xa are the signals received at station 1 and 2, respec- 
tively. Over the Integration period (typically 3 minutes) the model 
time delay Is approximated by 

A 

T - T + Tt (2,3-Z) 

where t Is a constant delay and T Is the delay rate. Maximizing the 
correlation function with respect to t and t results In the maximum 
likelihood estimates of delay and delay rate (Vlhltney, 197A]. The 
actual cross correlation process Is described In [Thomas, 1972a ,b; 
Whitney et al., 1976]. The statistical model for the estimation of the 
observables Is developed In [Whitney, 107A]. 

Precision estimates for the delay and the delay rate can be com- 
puted as a function of the system characteristics. However, they do not 
Include error sources such as the propagation mediunt, instrumental 
effects and modelling errors all described in section 2.4. The pre- 
cision for delay is given In [Counselman et al., 1979] as 

a(T) ^ (2.3-3) 

^^sp N 

Indicating that it is inversely proportional to the spanned bandwidth 

Af and the slgnal-to-nolse ratio (S/N) where 

sp 

/t t 

“ * a((j)) - 3.2 X 10^ ( — ^ I — i rpT ' seconds (s) (2.3-4) 

where 
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th 

the antenna diameter at the 1 atatlon (m) 

Af spanned bandwidth (Hz) 

sp 

number of tape-recorded and cross correlated samples 

th 

T the system temperature at the l'^" site (*K) 

®1 

Ej the antenna efficiency at Bite 1 

F the correlated flux density (Janskys)-(that fraction of 

the total flux density from the source that "survives" 
cross correlation) 

a(4i) Is the uncertainty In the estimation of the fringe phase, from 
each of the narrow separate bands. As an example for a typical Mark III 
3-mlnute observation, 

Di ■ D 2 ■ 50 m 

Ej a C 2 ■■ 0.5 

F "1 Jansky 

T - T - 100“K 

Si 82 

Af - 400 MHz 
sp 

Nj. “ 14(7.2 X 10®) bits (based on a 4 mblt/s sampling rate 
per track) 

we arrive at 0 (t) « 6.4 picoseconds (ps) . 

At the present state of the art, the delay precision ranges 
below the cm level for a 3 minute integration period. It can be seen 
from (2.3-4) that an Increase in the spanned bandwidth will allow 
deployment of a smaller antenna at one of the sites and not incur any 
loss in precision. The precision of delays Is Inversely proportional to 
the correlated flux density. However, on long baselines many of the 
compact sources arc partially resolved since angular resolution Improves 


with baseline length. This results In a decrease of F so that sources 
which show strong fringes on baselines of a few hundred kn become very 
weak on Intercontinental baselines. 

An expression for delay rate precision Is 

a(T) - s/s (2.3-5) 

Vn 

where t Is the total Integration time and Is the root-mean square of 
the sampling frequencies [Whitney, 1974]. 

2.4 Deviations from the Geometric Model 

The product of cross correlation Is a set of estimated delay 
and delay rates, and their precision estimates. The geometric observ- 
ables have been described In section 2.2.2. It Is left to describe 
those physical effects that cause the group delay and phase delay rate 
to differ from their corresponding geometric counterparts. These arise 
from Instrumental Imperfections, source structure, the propagation 
medium and other factors, all described briefly In this section. For 
more details, appropriate references are given. 

The frequency standards located at the various sites must have 
short- and long-term stability. The former Insures that the relative 
phase of the signals can be accurately recovered through cross correla- 
tion. The use of hydrogen masers effectively eliminates errors of this 
sort. The long-term stability of the clocks Is necessary In order to 
keep accurate time and prevent drifts In the relative clock behavior. 
This stability may falter at Intervals of time as short as eight hours. 
In this time period. If the long-term stability of the clock was 
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approximately 1 x 10" as can be achieved (or better) at present in 
the laboratory, this would lead to an error of 0.3 nanoseconds (ns) In 
time delay corresponding to an error of several cm in baseline length 
depending on the baseline chosen, if not corrected. Hydrogen masers, 
moreover, have been found in field work to be Influenced by atmospheric 
conditions and other environmental factors. Systematic patterns in the 
least-squares residuals may Indicate poor clock behavior. The usual 
remedy is to model these errors by polynomials as done in the next 
chapter. Other techniques include differencing of observations 
[Robertson, 1975] and the use of "clock stars" [Shapiro, 1979]. Anti- 
cipated technical Improvements in frequency standards and improved 
models will substantially reduce clock errors. See [Robertson, 1975] 
for a good example of how errant clock behavior is handled. The per- 
formance of hydrogen masers is discussed in [Vessot, 1979] and 
[Reinhardt et al., 1979]. 

Other Instrumental errors are caused by retardation and disper- 
sion of the signal as it passes through the cables and receiver compo- 
nents. These effects which can be of the order of several tenths of 
nanoseconds can be reduced significantly by phase calibration and cable 
measurement systems [Rogers, 1979]. 

Source structure Introduces unwanted noise (from the geodetic 
point of view) into the observables. The radio-sources are not gener- 
ally point-sources as assumed in section 2.2.2, and may exhibit compli- 
cated strvicture. Source structure maps are developed by radio astrono- 
mers which can be used to define a reference point for the source 
coordinates. Most of this information is derived by examining phase 
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closures around a triangle of stations since all other systematic 
errors cancel out. See [Hutton, 1976] and [Cotton, 1979] for more 
details on source structure. 

As In most geodetic systems, the propagation medium Is the ulti- 
mate limit on accuracy. The effects of the Ionosphere can be virtually 
eliminated by observing enough sources In two widely spaced frequency 
bands or by choosing a relatively high center frequency for which the 
Ionospheric effects would be small [Whitney et al., 1976]. These 
errors can be reduced to well under 0.03 ns In delay [Counselman, 

1976]. The dry component of the troposphere which Introduces an error 
In the time delay of up to 7 ns at the zenith can be modelled quite 
well based on recordings of surface metereologlcal data. In addition, 

It can be parameterized by a zenlttti distance thickness parameter scaled 
as a function of elevation angle [Ma, 1978]. The wet component of the 
troposphere poses the most serious problems though Its effect Is less 
than 1 ns In delay. The water vapor In the troposphere changes with 
respect to time and direction of observation. It Is hoped that with 
water vapor radiometry the total uncertainty In tropospheric error can 
be reduced from about 0.1 ns for the zenith direction to 0.03 ns. These 
errors map particularly into the vertical component of the baseline. 

As the accuracy of VLSI observations Increases and especially 
for longer baselines, relativistic effects must be considered. Electro- 
magnetic waves are deflected by the gravitational field of the sun 
according to Einstein's theory of general relativity, thereby affecting 
the time delay. For further details, see [Thomas, 1972], [Robertson, 
1975] and [Gourevltch et al . , 1979]. 
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Another effect Includes the gravitational flexure of large radio 
telescopes which changes the location of the VLBI antenna reference 
point [McGinnis at al., 1979]. For example, In the comparison of the 
Haystack-Westford baseline vector measured with VLBI and classical geo- 
detic methods there was a difference In the vertical component of 19 mm 
as compared to 2 and A mm In the two horlxontal components. 

By correcting for the gravitational flexure of the Haystack antenna the 
discrepancy In the vertical component was reduced to 6 millimeters 
[Carter, In press]. 

Inadequate geophysical modelling will also Introduce systematic 
errors Into the estimation process. These Include errors In nutation, 
precession, UTl and polar motion as well as Incorrect earth tide and 
ocean loading models. These effects will be discussed In more detail In 
Chapter 3. 

The adequate modelling or elimination of systematic effects will 
determine the attainable accuracies for geodetic and related parameters. 
This Is especially crucial for the detection of geodynamlc phenomena. 
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3. MATHmilCAL KODELS 


3.1 Introduction 

In this chapter, the various natheioatlcal models used In the 
VLBI Interactive Program (VIP) are described. In section 3.2, the mathe- 
matical models for the VLBI observables are derived. In section 3.3, 
singularity problems due to coordinate system definition, observ- 
ability conditions and critical configurations are summarized. Finally, 
In section 3.4, the radio-source observability equations are given. 

The choice of a parameter set for VIP was Influenced by optimi- 
zation studies related to the Polaris network. Therefore, the stress Is 
on earth orientation variation parameters. Of course, baseline para- 
meters are also of primary Interest. Source coordinates are needed in 
order to develop a reasonably accurate catalogue from which more accu- 
rate geodetic parameter estimation will follow. Clock parameters, 
though of no direct interest here, are necessary to make the ana3.ysls 
more realistic. Atmosphere parameters, though not Included In VIP, may 
be useful If metereologlcal data Is not sufficient [Ma, 1978]. Smaller 
effects that require long observational campaigns (e.g., geodynamlc 
phenomena, precession, nutation) have not been parameterized. 

The adjustment philosophy has been to avoid weighted parameters, rather 
to define estimable parameters which Implicitly supply the minimal con- 
straints needed for Invertiblllty of the normal matrix. All parameters 
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are eatlnated from the observations themselves without resorting to 
external information. 

3.2 Least Squares Adjustment Mathematical Models 

3.2.1 Introductio n 

In section 3.2.2 the "lne*’tlal" and terrestlal coordinate sys- 
tems are defined. The mathematical models for delay and delay rate 
observations are presented in section 3.2.3 and 3. 2. A, respectively. 

For each observable, the estimable parameters are defined. Section 

3. 2.5 is a description of the least squares algorithm. In section 

3.2.6 a simple model, suitable for covariance analyses, is presented 
for computing the correlation between delays observed simultaneously at 
a given epoch, from a multistation configuration. Possible model 
refinements are discussed in section 3.2.7. 

3.2.2 Coordinate Systems Definitio n 

In analyzing VLBI observations an "inertial" and terrestlal 
coordinate system need to be defined. In practice, a "nearly" 
inertial frame is defined with its origin at the s«1ar system bary- 
center. Tlie first axis is directed towards the mean vernal equinox at 
some reference epoch, conventionally 1950.0 and the third axis is per- 
pendicular to the mean equator and positive northward. The second axis 
completes a right-handed Cartesian coordinate system. The theoretical 
calculation of delay and delay rates are performed according to the 
laws of general relativity in this coordinate system [Counselman, 1976]. 
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Exprssslons for these observables are derived relatlvlstically by 
[Robertson, 1975]. Since arrival tines are measured by atomic 'docks 
at the various stations, they must be transformed to coordinate tine of 
solar-system bary centric coordinates [Robertson, 1975, appendix B]. 

The transformations from the geocentric origin to the solar-system 
barycenter is done using a planetary ephemerls. It should be noted 
that the use of the above coordinate system implicitly includes the 
effects of annual and diurnal aberration [Ma, 1978]. The reason for 
this coordinate system definition is to be able to easily combine VLBI 
observations with spacecraft tracking and Interplanetary radar data. 

In VIP, it is assumed that the source positions have been 
updated to their true-of-date coordinates at the initial epoch of 
observation (precession and nutation corrections are not applied in the 
program). In addition, it is assumed that the observables have been 
corrected for aberration and for relativistic effects. Therefore, the 
"inertial" coordinate frame is taken as a true-of-date geocentric sys- 
tem defined at the initial epoch of observation. 

The terrestlal (earth-fixed) coordinate system is defined with 
the X-axis directed towards the Greenwich mean astronomic meridian 
determined by Che Bill. The Z-axls is towards the average north terres- 
ital pole (the CIO pole). Tlie Y-uxls completes u right-handed Cartesian 
coordinate system. The origin of this system is arbitrary since the 
mathematical models only contain baseline coordinate differences. In 
the VIP experiments the station coordinates are taken in NASA's Space- 
craft Tracking and Data Network System (STDN) system. In practice, the 
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origin Is usually dsflnad by Chs sdoptsd coordinatss of ons VIBI 
antenna, given in sons terras tlal system. 

The reference orientation of the baseline vector with respect to 
the true*-of~date system must be defined externally at the initial epoch 
since VLBI observations are only sensitive to the relative orientation 
of the baseline vector as will be discussed in the next section. 


3.2,3 Time Delay Model 

The geometric delay was defined by (2.2-7) as 


‘*ijk " ’ ®k 

th 

which represents the inner product of the 1*^ baseline vector in the 
terrestial frame and the k*"*^ source unit vector transformed from the 
true-of~date system into the terrestial frame at epoch tj . Remember 
that the delay is given in units of distance. Adding a two term poly- 
nomial, whose coefficients Ac . and Ac . correspond to a relative off- 

01 ll 

set and rate, respectively, between the two clocks at the ends of the 
i " baseline, the delay can be written as. 


ijk 

where 6. 




+ + (3.2-1) 

is the Greenwich Apparent Sidereal Time (CAST) at 
epoch tj 

are the components of polar motion tiiat relate the true 

celestial pole ("Instantaneous" rotation axis of the 

earth) to the average terrestial pole at epoch J ((i^ is 

defined as positive along the Greenwich meridian and 

along the 270“ E meridian) 
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c th« sp««d of light 

tf th« Initial apoch of observation (in VIP taken as 0^ UT 
of initial day of observations) 

Th. Rj utrlc« r«pr...nt (rlghf h.nd) rotation, about th. aubacrlptad 
i axis by the angular arguoient in parentheses [Mueller* 1969]. The 
CAST, 0^ can be revritten as follow. 

0j - e« + w^uTij 

• 00 + W^[TAI-(TAI-UTC) - (UTC-UTl)]j (3.2-2) 

+ Eq. E. 


where Bp 

Eq. E. 
TAI 
UTC 
UTl 


CAST at 0 UT of the Initial day of observations 
equation of the equinoxes 
international atomic time 
coordinated universal time 

observed universal time corrected for polar motion 

W. conversion factor from universal to sidereal time, 

a 

In practice* UTC-UTl is interpolated from BIH Circular D five day 
values. For purposes of brevity, let us denote 


(UTC-UTl) 


J 


l*h 

at the j " epoch. Since and r|j are small quantities* expression 
(3.2-1) may be rewritten as 




1 0 

0 1 -n 

+ c[Ac^^ + ACj^(t^-t^)l 


cos6 

-aln0 

0 


sin0. 


CO80, 


cosdj^cosaj^ 

cos6j^sinaj^ 

sinfi. 


(3.2-3) 
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t 


H 

wh«r« AXj^,AYj|^,AZj^ «r« th« coordinate dlfferencea of the 1 baeeXlne 

In the terreatial frame 

are the true right aacenaion and declination of 
♦* 1 % 

the Ic ” aource, reapectlvcly. 

Equation (3.2-3) expreaaea the functional relationahip between 
the delay obaervatlona and the parametera Hated above. Of direct geo- 
detic interest are the baseline coordinate differences, AX. , AY., AZ, 

ill 

(from which the baseline length can also be determined) and the earth 
orientation parameters, hj* k ^ • Hie source coordinates, 
are of astrometric Interest. Eventually, their accurate dateiralnatlon 
will provide a catalogue of wall-distributed sources resulting in more 
accurate geodetic parameter estimation. The clock parameters, Ac^^^, 

ACj^ are nuisance parameters, defined in order to make the mathematical 
model more realistic. We will now examine which of the above para- 
meters are estimable. By estimabllity we mean that there exists a para- 
meter estimate which is unbiased, l.e., that the expected value of the 

/s 

parameter estimate should be equal to the parameter itself (E(X) « X) . 

In other words, the parameters can be estimated directly from the 
observables without introducing external Information (for example, para- 
meter weights). It follows that for an estimable parameter set (l.e., 
each parameter is estimable) the normal matrix (see below) is 
invertible. It is enough for one parameter to be non-estlmable for the 
normal matrix to be rank deficient (singular) , thereby preventing para- 
meter estimation. Using these properties, the estimable parameters cor- 
responding to the VIP mathematical models will be determined. 
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The normal matrix la by definition 


N - a’^PA 

where A la the dealgn matrix and P, the weight matrix of the obaerv- 
ablea . The elementa of A are the partial derlvatlvea of the obaervable 
with reapect to the correapondlng parametera of the mathematical model t 
In thla caae» the obaervable la the delay and the parametera of Intereat 
are 

AXf AY^, AZ^, Cj, rij, Ac^^, Ac^^ 

as described above. Equation (3.2-3) can be rewritten 

d^jj^ - -AX^[cos6j^cos(0j-Oj^) + 5jSln6j^] 

+Ayj^[cos6j^sln(e,-oij^) +tij8ln6j^] 

(3.2— A) 

-AZ^[sln6j^ - 5jcos6j^cos(0j-aj^) - njCos6j^sln(0j-Oj^) ] 

+c[Ac^j + Ac^j(tj-tj)] 

Taking the differential of with respect to the parameters listed 
above 

(3.2-5) 

the time delay 


■“Muk’ - ''ax + ''ay/"i *A2“"l 

+ + Aj^dd^ 

* ''Ac ''ao„'‘<''<'.1> 

OX IX 

where the A's are the required partial derivatives of 
with respect to the subscripted parameters as follows 
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( 3 . 2 - 6 ) 


" “CO86j^COB(0j-Oj^) - 


A ** 

^^i 

co86j^8in(0j-ct|^) + njSln5|^ 

(3.2-7) 

A ** * 

-[8ln6j^ - Cjcoa6j^cos(6j-aj^) - TijCosfij^slnCOj-Oj^) 3 

(3.2-8) 

Am ** * 

-AXj^sinfij^ + AZ^cos6j^cos(0j-Oj^) 

(3.2-9) 

A ■ 

AYj^slnSj^ AZ^coa6j^8ln(0j-oij^) 

(3.2-10) 

A *» 

W^co86j^tAX^sln(0j-otj^) AYjCos(Gj-aj^) 

- AZ^Cj8ln(0j-oij^) + AZj^njCos(0j-a^^)] 

(3.2-11) 

\ “ 

•\/«d 

J 

(3.2-12) 

A.. " 

sln6j^[AX^co8(0j-Oj^) - AYjj^8in(6j-aj^) 



-AZj,C.jCos(0j-aj^) - AZ^nj8in(0j“Oj^)l 
-co86^tAZ^ + AX^C^ "AY^rij] 

(3.2-13) 

A. " 

c 

(3.2-14) 

li 

c(tj - 

(3.2-15) 


If there exist linear relationships between the coefficients listed 
above, the column rank of the design matrix will not be full and the 
normal matrix will consequently be singular — Implying that not all of 
the above parameters are estimable. Neglecting the terms containing 
Cj , hj and being negligibly small, the following linear relation- 
ships are evident among the partial derivatives 



AX*A»„ *" 

1 i AXj^ 

(3.2-16) 


(3.2-17) 

+ AXjA^^ 1 

(3.2-18) 
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\ ' <3-2-w) 

Therefore, It is not possible to estimate all of the parameters of 
Interest from VLBI delay observations. In order to circumvent these 
rank deficiencies, a set of estimable parameters, closely related to 
the set listed above is defined which will allow the normal matrix to 
be inverted without the use of external information. 

Before defining this new parameter set, it is useful to present 
the geometric interpretations of the rank deficiencies, as expressed 
analytically by equations (3.2-16) to (3.2-19). The first three equa- 
tions show a linear dependence between various combinations of AX^^, 

AY^, AZj^ and h j » • These Indicate a rank deficiency of three due 

to lack of absolute orientation of the baseline with respect to the 
true-of-date frame which cannot be sensed by the observables. The 
origin of the terrestial system is arbitrary since the mathematical 
model is expressed in terms of coordinate differences. The scale, 
defined implicitly by the adopted speed of light, is Inherent in the 
observables. It is left to account for the rank deficiency expressed by 
(3.2-19). This is due to a lack of reference direction (origin of right 
ascension) for the true-of-date frame — the observables are insensitive 
to any absolute direction in inertial space. Thus, it can be seen that 
of the initial 10 parameters of interest only six may be estimated 
simultaneously (see 3.3.2). Notice that the clock of f set and rate para- 
meters are differences and not absolute. Therefore , any common errors in 
the epoch setting of the station clocks will be indlstingulshabxe from 
corresponding variations in earth rotation [Shapiro, 1979]. 
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Let U8 then define an estimable parameter set related to the 
original set. The earth orientation parameters will be redefined as 
follows. The total interval of observations is divided into several 
adjacent periods to be referred to as earth orientation steps (or steps) 
[Dermanls, 1978]. The three earth orientation parameters Dj* Kj 
will be rewritten as 


5 ^ - K, + 

- n 1 + (rij^-n |) (a>i) 


(3.2-20) 


where A refers to the step. The reference orientation of the base- 
line is defined by three parameters 5j, rii» <i referring to the average 
values of polar motion and UTC-UTlj, respectively, over the first step. 
For each subsequent step, a set of three earth orientation parameters 


^^lA “ ^A " 


AKja “ <A ” *^1 


(3.2-21) 


are estimated. They are Interpreted as variations in earth orienta- 
tion relative to the absolute orientation (implicitly provided by the 

tzh 

first step) averaged over the interval of time encompassed by the A 
step. These are the estimable earth orientation parameters and their 
estimates are influenced by the interval of time spanned by the first 
step and the number and spread of observations. By not including 
Hi, Ki in the parameter set, the linear relationships expressed in 
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(3.2-16)-(3.2-l8) have been broken without resorting to external infor- 
mation. This eliminates 3 of the 4 normal matrix rank deficiencies. 

The earth orientation variations can be added to rii» determined 
from other sources, for example, BIH Circular D interpolated values. 

For the purposes of VIP we can assume that 

Cl " Hi ■ - 0 , 

aJtliougli other values may be asf/igned in the program. 

A similar formulation will circumvent the fourth rank defi- 
ciency. The right ascension of one source will be constrained implicitly 
to its initial value by not Including it in the parameter set. We can 
write 

Oj^ - (k>l) (3.2-22) 

where ttj is the fixed true right ascension. This value will provide 
the reference, orientation of the origin of right ascensions. The cor- 
responding estimable parameters are the right ascension differences 
given by - Oj. Source right ascensions are non-estimable parameters. 
The declination of the reference source should be nearly equatorial to 
provide a strong definition for the reference direction. This can be 
seen by an examination of (3.2-12) since the right ascension partial is 

a function of cos6, . 

k 

This new set of estimable parameters is free of the rank defi- 
ciency of four exhibited by the initial set. Although the normal 
matrix is no longer singular, the estimation of baseline components is 
biased by any errors in the four parameters of orientation a^, T)j, 

Kj as will be shown below. From this point of view, the baseline 
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components AX^> AY^. AZ^ are non~estlmable parameters and again we shall 
resort to defining a corresponding set of estimable ones, 

[Arnold, 1974] • respectively, according to the following derivation. Let 
us rewrite (3. 2-S) , using (3.2-20) , (3.2-21) and (3.2-22) In terms of the 
estimable parameters discussed above, neglecting terms containing 0 
and K 


+ A„ tdAY. -AX.da. d-AZ.dn. +H.AX.dK,l 
i 11 II al* 

+ [dAZ^ + AX^d5j - AY^dm] 

1 


(3.2-23) 
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where the partial derivatives (the A’s) correspond directly to those 
given In (3.2-6) to (3.2-15). The partial derivatives of T^, e^, 
correspond to those of AX^, AY^, AZ^, respectively. The differential 
relationships between these tv'o sets are given by the bracketed terms 
in (3.2-23) 


dX^ “ dAX^ - AZ^dC^ + AY^dgj 
de^ » dAY^ -f AZjdtij - AX^dpj 
do^ - dAZ^ + AX^dCj - AY^dPj 


(3.2-24) 


where 


(3.2-25) 


d0| - da, - WjdK, 


Implying that these two rotations are Inseparable. The differential 
relationships between the parameters can be re-written in matrix form 
as 


, 

o. 

H 

H- 

1 


dAX^“ 


— 1 
1 

de^ 

- 

dAY^ 

+ R^(d5^)R, (dnj)R,(d3j) 

AY^ 



_dAZ^_ 


AZ^ 

m m 



'dAX^' 


0 

de, 

-d5,‘ 


1 — 
■H 
X 
< 

1 

- 

dAYi 


-de^ 

0 

dn, 


AY^ 


dAZ^ 



-dn, 

0 


AZf 


(3.2-26) 


where d?,, dn^, dg^ are errors in the initial reference orlenatlon 
assumed to ba of small magnitude. are the rotation matrlce. 
described earlier. Of course, the smaller these errors the more 
closely T^, e^, will "resemble” the baseline components. The 
Importance of accurate initial orientation parameters is especially 
apparent for long baselines. For example, from (3.2-24), for a base- 
line with ■ 4000 km, an error d?, « OVOOl will contribute to a 
change of 2 cm in the "estimated" AZ component (see Appendix B.l). 

Baseline lengths, on the other hand, are estimable quantities 
being unbiased by the errors in the reference orientation. This can be 
shown by writing the baseline length, as 

^1 " (3.2-27) 
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Then, 




AX^dAX^ + AY^dAY^ + AZ^dAZ^ 


(3.2-28) 


Substituting (3.2-24) into (3,2-28) yields 

dA. 

- AX^(dT^-AY^da,+AZj^dCi+ W^AY^dKj) 

+ AY^(de^ + AX^da,~AZj^dn,- W^AX^dic,) 
+ AZ^(d 0^ - AX^dC , + AY^dn , ) 


thus, 


^-AX^dT.+ iVjde. + AZjd^ 


(3.2-29) 


(3.2-30) 


Comparing (3.2-28) and (3.2-30), it follows that is unaffected by 
errors in ai, Ci, Hi and <i which is obvious since distance is invari- 
ant of coordinate system definition. However, baseline lengths as well 
as components will vary due to earth tides and geodynamlc phenomena, and 
therefore these phenomena may be parameterized as will be discussed in 
section 3.2.7. 

In VIP, the baseline length standard deviations are estimated by 

propagation of errors from the baseline "components" t . , e . , 0 . . The 

ill 

mathematical model is given by equation (3.2-27). The 
variance-covariance matrix for distances, is given, using the nota- 
tion by [Uotila, 1967] as 






(3.2-31) 
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where G Is the matrix of partial derivative* of with reapect to each 
component. covariance matrix of the baeellne 

"components" retrieved from their corresponding elements In the 
variance-covariance matrix of estimated parameters. 

It Is appropriate to summarize the previous discussion by 
listing the estimable parameters recoverable from delay observations 




k 


baseline "components" contaminated by errors In the 

reference orientation 

baseline distances 

source declinations 
right ascension differences 


polar motion variation components 


UTl-UTC variations 

ACjij^.ACjj^ relative clock offset and rate, respectively. 

3.2.4 Time Delay Rate Model 

The geometric delay rate was defined In section 2.2 as the time 
derivative of the geometric delay. Including the clock parameters the 


delay rate is modelled 




• \ * “‘.1 


(3.2-32) 


Differentiating (3.2-4) with respect to time 


d^jj^ - u^coa6j^ {AX^sln(6j-aj^) + AY^co8(0j-«j^) 
- AZj^t5sln(ej-aj^) - ncos(ej-aj^)]} 


(3.2-33) 


+■ cAc, 




where 


e dt 


is the spin rate of the earth, S) thr Inetantaneoue earth rotation vec- 
tor. The magnitude of the terms containing ^ and n In (3.2-33) are 
negligible, Indicating that the delay rate is effectively Insensi- 
tive to the AZ component of the baseline. It follows that only the 
length of the equatorial projection of the baseline can be eotlmated. 
In addition, the delay rate Is unaffected by clock offset variations, 

Ac .. Furthermore, examining (3.2-32) 

01 


X B 


is orthogonal to n and, thus, the origin of declination is undefined 
[Counselman, 1976] as well as the right ascension origin. The discus- 
sion of the parameters estimable from delay rate is identical to that of 
delays except that in this case AZj^(o^) and Ac^j^ are deleted, and decli- 
nation differences replace 6j^. Thus, an expression similar to 

(3.2-23), corresponding to delay rates 


d(di^kil> “ ^*^^1 " ^^1*^^! ^^1*^^ 1 ^ 


+ Ag [dAY^ + AZj^dn , - AX^d0 j ] 

(+ [dAZ^ + AX^dSj-AY^dnJ)’ 
i 


(3.2-34) 




*Negligible. 
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wh«r« 6 , l 0 th« daclimitlon implicitly constrained to Its a priori 
value by not Including It in the parameter set. All the other terms 
have been defined in section 3.2.3. The partial derivatives of the 
delay rate with respect to the subscripted parameters are 




(3.2-35) 

(3.2-36) 


(A 


0 ^ " 




(3.2-37) 


- ^iIC^co8(ej^-a^) + njj8in(0jj^-O,|^)]} 

^V?,) " (3.2-38) 

(3.2-39) 

(3.2-40) 


^^i^Vin(0jr\) “nj^cos(0 -oj^)]} 


(3.2-41) 


Ac 


li 


(3.2-42) 


From (3.2-41) it is evident that the delay rate la Ineenaltlve to the 
declinations of sources near the equator. 

The delay rates are less Important than the delays because of 
their relatively lower accuracy and reduced estimable parameter set. 

* 


Negligible - not Included in VIP as well as all other terms 
JariJfS" “<■ .l»Uarly fo^™: delay 
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How«var, delay rate obeervatlona do have the advantage of being 
unambiguously eatimated and, thus, may be estimated with rela- 
tively simple equipment. In addition, Fanselow [1978] states that 
the delay rates aid in reducing correlations between certain 
parameters . 


3.2.5 Adjustment Algorithm 

The adjustment algorithm used in VIP is the standard method of 
observation equations of the form [Uotila, 1967] 

L - P(X ) 
a a' 

where L is the theoretical value of the "observed” quantities, delay 
and delay rate, related functionally to the theoretical values of the 
parameters. The function F is given by equations (3.2-3) and (3,2-33) 
for delay and delay rate, respectively. The non-linear function F, in 
each case, is linearized by retaining the first-order term of the Taylor 
series expansion about the approximate values of the parameters, Xg, 
such that 

3F(X ) 

- F(X.) + — (X_^-X,) 

- Lo + A X 

where Lg ■ F(Xg) is the vector of approximate values of the observed 
quantities based on the approximate parameter vector, Xg and computed 
from equations (3.2-3) and (3.2-33). The design matrix of partial 
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3F(X ) 

derivatives A " — Includes the elements given by equations 

* X * X, 
a • 

(3.2-6) - (3.2-15) and (3.2-35) - (3.2-42). X - X - X Is the vector 

a 0 

of parameter corrections to be applied to the approximate parameter 
estimates, X , to yield X • the adjusted parameters. The theoretical 

0 A 

observable, can be separated Into the actually observed quantity vec- 
tor, (In this case group delay and phase delay rate estimated from 
the cross correlation process) and the vector of residuals, V, resulting 
from observational errors. Then, 

Lj^ + V-Lj+AX 

V - A X + L 

where L ■ 

T 

By minimizing the sum of the squares, V PV, the least squares estimate 
for the parameter correction vector, X is 

X - -(a^pa)“^a’^pl - -n"^u 

where P Is the Inverse of the variance-covariance matrix for the observ- 

2 

scaled by 00 , the a priori variance of unit weight. 

The a priori covariance matrix of the parameters Is given by, 

- 05 (A^PA)"^ 
a 

The matrix is the basis of the VIP covariance analysis. The 
a 

a posteriori covariance matrix Is given by 


ables, 2 

h 
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whfltrc 


A 

z. 


'*'2 T -1 

a;(A^PA) ^ 


V^PV 
n - u 


''2 

Og la the a posteriori variance of unit weight, n la the number of 

T 

observations and u, the number of parameters* The scalar V PV can be 
computed from 


T T T 

V^PV - L^PL + XU 


therefore, there Is no need to compute each residual. However, In 
practice the residuals usually contain Information on systematic 
effects, especially errant clock behavior. Since VIP Is mainly 
Intended as a covariance analysis program, the residuals are not com~ 
puted when the least squares solution option is specified. 

The VIP least squares algorithm uses the equations listed 
above. The normal matrix, N Is filled In a sequential manner and in 
upper triangular form in order to conserve on storage requirements. 

This is crucial on TSO where the limit Is 256K. Triangular storage 

2 

requires u(u+l)/2 storage locations as opposed to u In the full case. 
No attempt Is made to exploit normal matrix sparsity patterns although 
this may become necessary for larger parameter sets. VIP Is dimen- 
sioned to accept a parameter set of size 62 although this could be 
Increased up to the storage limit of 256K. In order to simplify dimen- 
sioning all matriceii ire stored in vector form using the SSP subrou- 
tine LOG for bookkeeping purposes [IBM, 1970]. The SSP routine, DSINV 
which handles matrices stored in upper triangular form is called to 


invert the normal matrix. 
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since simultaneous observations from several stations to a given 

source at a particular epoch J are correlated (as described In the next 

T 

section), the N and U matrices and part of V PV are filled epoch by 
epoch as follows 


U - E A/P.L. 
j-1 J J 

V^PV " X^U + E l/p.L 
J-1 J J J 

where E Is the number of simultaneous observation sets, each set con- 
taining the observations of one epoch. The P^ portion of the weight 
matrix is block diagonal, each block having Its dimension equal to the 
number of Independent baselines observing simultaneously at that epoch. 
This will become clear in the next section. The above summations 
assume that observations at different epochs are uncorrelated which Is 
in accordance with the VIP mathematical model. In practice, such obser- 
vations may be correlated but only as a result of unmodelled systematic 
effects such as those resulting from the propagation medium. 

After Inversion of the normal matrix, and multiplication by the 
variance of unit weight, the estimated standard deviations of the para- 
meters are computed by taking the square root of the diagonal elements 
of the resulting variance-covariance matrix of parameters. In addition, 
the correlation matrix of parameters is computed from 
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'x.x 


Ov cr„ 
^1 


where o„ „ Is the covariance of parameters X. and X. and a„ and 

Xj^Xj i j 

are their respective standard deviations. The correlation matrix 
describes the Interrelationships among the parameters. A value of 
|p„ „ 1 close to unity Indicates that the parameters are highly depend- 

XiXj 

ent while a value of unity Indicates a singularity and complete linear 
dependence. High correlations may result In Ill-conditioned matrices 
and thus unstable systems whose solutions are circumspect. 

Ill conditioning of the normal matrix Is reflected by the ratio 
of the largest and smallest eigenvalues. They are competed in VIP 
using the SSP routine, DEIGEN, which outputs the eigenvalues In 
descending order of magnitude. A relatively large ratio will Indicate 
111-condltlonlng possibly resulting from a critical geometric configura- 
tion (see Section 3.3.3) « 



3,2.6 Weighting of Observations 

VLSI observations are usually performed simultaneously from all 
participating stations unless mutual source visibility makes this Impos- 
sible. In accordance with the VIP mathematical model, simultaneous 
observations to a particular source at a given epoch are correlated. A 
simple model, suitable for covariance analyses, for computing these cor- 
relations will be described below. This formulation assumes that the 
delays are all observed with equal precision, a reasonable assumption 
for covariance analyses. Typical precisions are 0.1 us (3 cm) for 
delay and 0.1 ps/s (0.108 m/hr) for delay rate. 
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The following discussion will address a triangle of stations for 
the sake of description but can be extended to any closed configuration. 
As described In Chapter 2, the raw observables are the bits recorded on 
magnetic tapes at the three sites. The d*.x4y (and delay rate) Is esti- 
mated by cross-correlation of the tapes. Denoting the time delay 
between stations 1, J as It follows from the mathematical model 
that, 

t, 2 + Ta, + T,| - 0 . (3.2-43) 

Thus, any one of the delays Is linearly dependent on the other two. In 
other words. If two delays have been estimated then, theoretically, the 
third one Is completely determined (Shapiro, private communication) and 
does not provide Independent Information. In this example, there are 
three possible combinations of two Independent delays. Regardless of 
the chosen combination, the parameter estimates should be Identical 
since all three sets of tapes, containing the same Information In any 
case, are required. If the correlations between the observables, at 
each epoch, are neglected, there will be three different sets of para- 
meter estimates, one for each condblnatlon. 

The delay, conceptually. Is the difference In times of arrival 
of a given portion of a wavefront at two antennas. Therefore, In tri- 
angle 1-2-3 the delays for one epoch can be written as 

‘’f 12 ■ ^2 - 1 1 

Tjs “ t3 - t2 (3.2—44) 

'^31 - 
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In matrix form 


“tiz" 


1 

0 

1 ! 


■ 

T 2 I 

■1 

0-11 



-Til. 


1 

f— 4 
1 

0 

r-i 

1 


-^1 - 


G T 


(3.2-A5) 


Assume that the variance-covariance matrix of the "observed" times 
Is a diagonal, l.e., that all observations are of equal precision. Let 
us further assume that It Is the Identity matrix since at this point we 
are Interested solely In the correlations between delays. By propaga- 
tion of errors, the variance-covariance matrix of observed delays for 
one epoch of observation Is 


Z 


T 


T T 
- GG^ 

■ 2 -1 - 1 ‘ 

-1 2 -1 

.-1 -1 2 . 


(3.2-46) 


However, the determinant of this matrix (of rank 2) is zero and thus cannot 
be Inverted. This Is just a restatement In mathematical terms of the fact 
that the three delays are dependent . Clearly, parameter estimation Is 
impossible In this case and one delay must be eliminated. It makes no dif- 
ference which one since, using this model, the parameter estimates and 
their variances will be identical using any two of the time delays. Let 
us choose Tj 2 and T 23 . Then, for one epoch of simultaneous observations 



49 


disregarding the scale factor for the moment. Scaling this matrix to 3 cm 
precision In units of distance (equivalent to 0.1 nanosecond) and thus 
replacing x by d 

0.0009 -0.000A5 

® -0.000A5 0.0009 

This matrix Is then inverted and the first element (upper left-hand 

corner) Is scaled to unity, the scaling factor being the a priori 

2 

variance of unit weight, 0o* The complete weight matrix in this 
example is 2 x 2 block diagonal 


1 


1 

2 


P 








(3.2-47) 


2 

where Oo ■ 0.000675. As with all other arrays, the weight matrix is 
stored In upper triangular vector form. 

Correlations between simultaneously observed delay rates are 
computed in the same manner. In analyses Involving both observables it 
is assumed that delay and delay rates are uncorrelated. As mentioned 
earlier delays at different epochs are assumed to be uncorrelated, and 
similarly for delay rates. 

The above discussion indicates the importance of including the 
proper correlations among simultaneous observations. Otherwise, para- 
meter estimation is not unique and it is meaningless to perform a 
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covariance analysis. Using a diagonal weight matrix (neglecting the 
correlations) with all three time delays will yield a unique set but 
with parameter precision estimates that are overly optimistic. This 
may not be significant In the triangular configuration described 
above. In general, however, for an N-statlon configuration there are 
N(N-l)/2 possible baselines (tape combinations) but only N - 1 Inde- 
pendent ones. For example, in a six station configuration there are 
15 possible baselines, only 5 being Independent. 

The weighting procedure described above Is highly simplified 
but appropriate for covariance analyses. In analyzing real data, the 
matrix is much more difficult to determine and the weight matrix Is 
taken as diagonal. However, unless the true correlations are known, 
the least squares estimates may be quite misleading especially in 
larger networks. 

3.2,7 Model Refinements 

The mathematical models described in sections 3.2.3 and 3.2.4 
are suitable for the type of applications for which VIP is Intended. 
The parameter set chosen for VIP was Influenced by studies of the 
Polaris triangle, i.e., monitoring of earth orientation variations. 
Other effects such as nutation, precession, crustal movements, earth 
tides and ocean loading were not included. In the handling of real 
data, though, these other phenomena must either be parameterized or 
compensated for by a priori information in order to correct for their 
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Influence on the observables. Otherwise, the estimated parameters would 
be contaminated by their effects. 

Robertson [1975] has estimated the precession constant, the rate 
at which the Earth's spin axis rotates about the ecliptic pole, from 
VLB! observations spread over approximately four years. In simulation 
studies, Dermanls [197''] defined three rotation angles to model the 
total effects of precession and nutation. A step approach similar to 
that described In section 3.2.3 for earth orientation was used since 
only relative variations may be sensed by the observables. 

Robertson [1975] was able to estimate the Love number, h, 
related to radial displacements caused by the tidal potential. Since a 
time delay can be estimated every few minutes, this provides Ideal con- 
ditions from the point of view of earth tide analysis [Bonatz et al., 
1978]. 

Geodynamlc phenomena may be estimated from VLBI observations, 
by observing relative changes In the baseline components. When Involved 
In a long observational campaign many data sets are generated. Although 
baseline components are non-estlmable, the adoption of one reference 
baseline orientation for all the data sets at least will Insure that the 
estimated x, E, a parameters will refer to a consistent coordinate sys- 
tem. In this case, errors In reference orientation will cancel out. 
Otherwise, the differences In these parameters due to the varying 
reference orientation will look like tlme-like variations, although in 
reality they will only be due to inconsistent coordinate system defini- 
tion. 
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3.3 Singularity Problema 


3.3.1 Introduction 

In the least squares process, when the normal matrix has rank 
less than Its dimension It Is singular and cannot be Inverted, thereby 
preventing parameter estimation. In the VLBI case, more specifically 
using the models of VIP, singularity problems may occur for a variety of 
reasons. In this section we will review these problems. First, It 
should be noted that by defining the estimable parameters In 3.2.3 and 
3.2.4, singularities due to coordinate system definition have been elim- 
inated. As discussed previously, the origin of the terrestlal system Is 
arbitrary and the scale Is Inherent In the observations themselves. The 
reference orientation of the terrestlal frame with respect to the 
true-of-data frame must be specified and this is done by parameterizing 
the earth orientation parameters as variations relative to the values 
assumed for the first step. In addition, the singularity due to lack of 
a reference direction for the true-of-date frame Is eliminated by esti- 
mating right ascension differences relative to one fixed right ascen- 
sion. For delay rate observations only, the origin of declinations must 
also be specified. 

3.3.2 Observation Singularities 

In least squares estimation the number of observations must, of 
course, exceed the number of parameters. In VLBI, these observations 
must be distributed correctly over a minimum of three sources, other- 
wise a singularity will occur. This can be seen from the following 
analysis that has been performed previously by Robertson [1975] and 
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described very clearly by Shapiro [1978], both for a smaller parameter 
set than Included In VIP. 

Equation (3v2»4) can be rewritten as (dropping the subscripts) 

d - K,cos(0-a) + K 2 sln(e-a) + K, + K.,t (3.3-1) 

where Kj ■ (-AX + AZ^)cos6 
Ka - (AY + AZn)cos6 
Ks - (-AZ - AX? + AYn)sln6 + cAc© 

Ki| ■ cAci 

For a given baseline the K terms are constants to a first approxima- 
tion, but vary slowly with respect to time due to polar motion varia- 
tions (as well as precession, nutation and earth tides) . The terms 
Kicos(6-a) and Kzsln(6-a) are both diurnal sinusoids, remembering 
though that 0 Is affected by UTl-UTC variations (the K term). The 
amplitude of these sinusoids given by K| and K 2 are functions 
of the baseline vector, and the source declination. The two 
curves are shifted In phase by 90“ since sln(6-a) ■ cos(0-a--2)* Tli® 
angular frequency of the sinusoids Is given by the rotation rate of the 
earth. This can be seen by expressing 


where (ft le the phase of the sinusoids relative to some Initial epoch, 
is the rotation rate of the earth, and 6 Is the Greenwich sidereal 
time. The sum of these two sinusoids Is again a sinusoid of the general 
form 

Kcoa(^-4<i» t) 

¥ 0 

idiere Is the resulting phase and K, the amplitude. Therefore, 

d - Kcos((j».4tc t) + K, + K^t (3.3-3) 

which represents a straight line added to a diurnal sinusoid. 

Assume that delay observations are performed from one baseline 
to one source. From the discussion of section 3.2.3, over the first 
step the following parameters are estimable: Tj, e^, 6,, Ac^j, 

Acjj^, a total of 6 parameters. An examination of (3.3-3) Indicates, 
though, that only 4 Independent parameters K, 

be estimated from at least 4 observations to one source. Three addi- 
tional observations to a second source will enable 3 more Independent 
parameters to be estimated, another set of K, (j>g, Kj~a total of 7. 

Note that is common to observations of all sources. Two new para- 
meters, Oj”®! ^2 added to the set of interest— a total of 

8 parameters. Thus, observations to two sources still yields a singu- 
lar case with respect to the parameters of Interest. In a similar 
manner, 3 additional observations to a third source will allow estima- 
tion of ten independent parameters. In this case Olj-Oj, 6j will be 
added to the set of Interest — a total of 10 parameters. Thus, over the 
first (earth orientation) step we are able to estimate 
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T|, C|i 0 |, ot 2 **(X|( Oj-cii, 6|» 6j» Acg|» Ac 1 1 

from at least 10 observations distributed as described above to 3 
sources. For each subsequent step, we can estimate the 3 earth orien- 
tation variations AC^j^, Ak^^ as described in 3.2.3. As shown in 

the next section, only two of three of these parameters are estimable 
from observations from one baseline. In this case for each extra step, 
two more observations per step will be required to any of the three 
sources observed over the first step. In multl-basellne configurations 
all 3 earth orientation parameters may be estimated. In addition for 
each extra baseline the parameter set Increases by five, e^, o^, 

Ac Ac .. Thus, observations to the minimum 3 sources must be 
increased accordingly Of course, for the sake of redundancy the number 
of observations always exceeds the minimum number required (increase In 
the degrees of freedom) . In addition, for Improvement of the geometric 
strength of the observations, more than the minimum three sources are 
observed . 

For a similar analysis of delay rate observations, see 
[Robertson, 1975]. It is important to note that the addition of these 
observations do not add any independent information whereby the number 
and distribution of observations could be reduced. As mentioned 
earlier, delay rates have a reduced parameter set associated with them, 
though adding redundant information but of relatively lower quality. 

3.3.3 Critical Configurations 

There remains one additional category wherein the normal matrix 
is rank deficient and this can be classified as critical 
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basallne configurations. Basaline orientation approaching these special 
cases will result In high correlations between certain parameters and 
111-condltlonlng of the normal matrix. 

It will be shown here that observations from a single baseline 
are sensitive to only two of the three earth orientation variation para- 
meters* AC, An, Ak. To understand this let us first examine the pos- 
sible critical configurations of the one baseline case. Consider a 
baseline parallel to the earth's axis of rotation observing a source at 
"infinity.” Examining the partial derivative (3.2-11) It is evident 

that since AX^ - AY^ - 0 (and neglecting the terms containing C and n) 

1 1 

A* “A “0 

Ak k 

and, therefore, che delay (and rate, see (3.2-37)) is Insensitive to the 
UTl-UTC parameters. Attempting to estimate these parameters will result 
in a singular normal matrix. Consider a baseline parallel to the equa- 
tor whose midpoint Is situated on the Greenwich meridian or at 180" 
longitude. In this case AX^ * AZ^ ■ 0 implying from (3.2-9) that 

so that for this configuration the delay (and rate, see (3.2-38)) is 
Insensitive to the A? parameters of polar motion variation. Similarly 
for a baseline parallel to the equator and whose midpoint is at 90"E 
or 270"E longitude from (3.2-10) 


A._ » A “ 0 

An n 
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and the delay (and rate, aee (3.2-39)) le ineeneltlve to the An com- 
ponent of polar motion. For example, continental United Stutea 
eaat-weat baseline observations are hi.rdly sensitive to the An para- 
meters. It is now evident that observations from one baseline can be 
used to estimate only two of three earth orientation parameters inde- 
pendently since a change in the orientation of one baseline is com- 
pletely described by two distinct rotations. The choice of which two 
to choose will be dictated by the orientation of the baseline and 
examination of the magnitude of the partial derivatives of (3.2-9), 
(3.2-10) and (3.2-11). Note that the addition of any number of parallel 
baselines to any of the critical configurations listed above will not 
eliminate the rank deficiency. However, observations on any two 
non-parallel baselines will allow estimation of all three earth orienta- 
tion variation parameters (over each step — except the first) since a 
change in orientation of a plane in space is fully described by three 
independent rotations. An exception to this is given in the next para- 
graph. 

A baseline parallel to the equatorial plane (observing delays) 
constitutes another case of a critical configuration. In this case, 

ASSj * 0 and from the partial derivatives (3,2-6), (3.2-7), and (3,2-13) 
the following linear relationship is evident 

(3.3-4) 

assuming the terms containing 5 and n are negllg'^ble. This wfll result 
in a rank-deficient normal matrix. In geometric terms, the origin of 
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declination la not seiiaad by the obaervatlons alnca the baaallna Is 
orthogonal to tha Instantaneous earth rotation vector, This Is 
similar to case of delay rates described In 3 2.4. Therefore, for such 
a configuration, the parameter set must be modified by Introducing 
declination differences as parameters Instead of declinations. It 
should be noted that for any number of east-west baselines the con- 
figuration will be critical when estimating the regular delay parameter 
set described In 3.2.3. In practice, baselines parallel to the equa- 
torial plane are not very common; however, observations from a baseline 
approaching this configuration may result In an all conditioned system 
and high correlations between certain parameters. The ratio of maximum 
to minimum normal matrix eigenvalues Is a good Indication of an 
ill-conditioned system. For typical "non-critlcal" VLBI baseline con- 

5 6 

figurations this ratio is of the magnitude 10 or 10 . For a 
near-critical configuration this ratio will be several orlers of magni- 
tude larger. 

3.4 Source Visibility Equation 

In planning an observation schedule the first factor to be con- 
sidered Is which sources are visible at a particular epoch and from 
which stations. This information is displayed in the visibility matrix 
computed for each source. The dimensions of each matrix are N x 24 
where N is the number of stations and the columns refer to the epoch of 
observation at one-hour intervals. The elements of the matrix ate out- 
put by VIP as zenith distances when the source is visible, a double 
asterisk when not. An example is given In Appendix B. 
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The source unit vector s in the inertial frame Is given aa 


s 


"cos 6 cosa' 
cos6 slna 
sir'' 


The station unit vector ^ Is given In the terrestlal frame as 


X 


’eo8(}) cosX"* 
cos4> slnX 
.sini() 


(}) and X , the geodetic latitude and longitude of the station, 
respectively. The zenith distance is given by 

z " cos ^(5*q) (3.4-1) 

where q corresponds to s rotated Into the terrestlal frame by Rs(6), 

0 being the Greenwich sidereal time (see eq. (3,2-3)). Polar motion 
has beet', neglected. Any cutoff angle may be specified for the accept- 
able zenith distance which is usually taken as 80" (or less) because 
of the large '.efractivity effects for observations near the horizon. 


4. CONCLUSIONS AND FUTURE RESEARCH 


A*1 Summary and Conclusions 

A VLSI covariance analyals Interactive Program (VIP) Is pre- 
sented In Appendix A as a tool for experiment planning, simulation 
studies and optimal design problems. Explanatory tables, figures and 
the necessary JCL are Included for ease of adaptation and operation. 

The sample session Included In Appendix B, consisting of two experi- 
ments, Illustrates some of the capabilities of the program and the 
advantage” of working In the Interactive mode. The program Itself Is 
well documented In case the user wishes to Incorporate his own modifica- 
tions (e.g., expanding the parameter set). By an explanation of the 
theory on which the program Is based and of the mathematical models which 
It Incorporates, an overview of the VLBI process Is given. 

The Introductory chapter touches upon the past, present and 
future aspects of VLBI as well as Its applications to geodesy and the 
related fields of astronomy, geophysics and geodynamics. 

Chapter 2 opens with a description of the basic VLBI geometry 
and of the quantities that are of primary Interest — the time delay and 
time delay rate. Their estimation requires expensive and sophisticated 
Instrumentation used for data collection and cross-correlation of the 
recorded tapes. The basic components of the system have been described 
(In broad terms) as well as their functions. The point Is made that to 
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arrive at geodetic parameter estimates, a two-stage estimation proce- 
dure Is required. The first one results In maximum likelihood estimates 
of group delay and phase delay rate and their corresponding precision 
estimates from cross-correlation of the recorded tapes. The basic 
observables In this stage are the raw bits recorded on the various 
tapes. The second adjustment, by least squares. Incorporates the above 
estimates as observables to estimate the relevant geodetic parameters. 

Of course. If the full covariance matrix of the "observables” (delay and 
delay rate) is available, there are no problems with this two-step 
procedure. However, In practice, this Is not the case (see Section 
4.3). In addition, unmodelled systematic effects further Influence the 
estimation process and tend to reduce the reliability of the final esti- 
mated parameters. Clearly more research is needed to improve the mathe- 
matical models and for a more rigorous statistical treatment of the 
data. 

In Chapter 3 the mathematical models used in VIP are described. 
The parameters estimable from delay measurements are derived. Except 
for baseline lengths and source declinations all estimable parameters 
are variational in nature, relative to the initial orientation of the 
inertial and terrestrial frames. These parameters Include baseline 
components that are contaminated by errors in the reference orientation, 
polar motion and UTl-UTC variations. The Importance of adhering to a 
standard reference orientation when combining several data sets is 
stressed so that all parameters will refer to a unique coordinate sys- 
tem. The delay rate model Includes a reduced parameter set. The third 
"component" of the baseline Is non-estlmable and only declination 
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differences may be estimated. Next, the VIP adjustment algorithm, 
the sequential observation equation method. Is described. A model 
Is presented for determining correlations between simultaneous obser- 
vations used to formulate the variance-covariance matrix of the 
observables. It Is suitable for covariance analyses and Illus- 
trates the importance of developing a model to handle real data. 

Next, possible model refinements are indicated Including parameteri- 
zation of precession and nutation, earth tides and geodynamlc phenomena. 
Finally, singularity problems associated with the models described In 
this report are summarized. These occur from incorrect distribution 
and number of source observations, and from critical baseline con- 
figurations. 

In the next two sections, two future areas of research will be 
discussed. Each one will make use of VIP and are actually the Impetus 
for Its development. 

4.2 Optimal Design Problems 

VIP is Intended as an aid in VLBI simulations by providing 
a priori lower bounds on the expected variances of baseline and earth 
orientation parameters to be estimated from a given experiment. These 
simulations may Include first- and second-order optimal design problems. 
Design problems. In general, may be approached from philosophically dif- 
ferent but related points of view. In experiment planning, a parameter 
(or set of parameters) Is required at a certain level of accuracy and 


the problem is to determine the conditions necessary for this require- 
ment to be met, if it can be met at all. In axperlment simulation, 
given a set of preliminary concltlons we wish to determine the upper 
bounds on accuracy for a given parameter set. VIP can be of help in 
studying both approaches although it is more directly amenable to 
experiment simulation. 

VIP was developed while studying the proposed Polaris triangle 
(Westford-Ft. Davls-Rlchmond) of NGS to be dedicated to the monitoring 
of earth orientation variations. Therefore, the description of 
first- and second-order design problems will use this network for 
explanatory purposes. First-order design may be defined, generally, as 
the selection of station sites for the ''optimal" estimation of a given 
geodetic parameter set. The criterlor, for optimality (in any sort of 
design) may vary, one example being the minimization of the trace (or 
partial trace) of the variance-covariance matrix of parameters. 
First-order design may be approached according to one of the viewpoints 
described above, for example: 

1) Suppose we wish to determine 24-hour averages of polar motion 
and earth rotation variations to the 10 cm and 1 ms level respectively. 
What are the minimal conditions necessary, e.g., number and choice of 
stations (given a source catalogue) for those accuracies to be met? Or 
alternatively, 

2) Suppose there are five available stations for the Polaris 
network but only three may be chosen. Which three would allow the 
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"beat" overall determinations of earth orientation variations? Obvi- 
ously, the first question Is more difficult because of Its absolute 
nature, while the second one requires only a relative answer. The 
presence of unmodelled systematic effects In the actual measurement 
process (particularly at the present VLSI state-of-the-art) may make 
good a priori accuracy estimates (via covariance analyses) difficult to 
obtain, especially for question 1 above. 

Ma [1978] studied the first-order Polaris network design prob- 
lem from an experiment simulation point of view. He Introduced the 
effects of typical systematic errors by including model error para- 
meters in the covariance analysis. The first-order design problem has 
not been of primary Importance since antenna availability in conjunction 
with economic, and political considerations have almost totally con- 
strained its solution. However, with the development of portable 
antennas and allocation of greater resources, this problem assumes 
greater relevance. Dermanls [1977] derived parameter sensitivity vec- 
tors in studying first-order design for earth orientation and baseline 
parameter estimation. 

The second-order design problem may be defined as follows; 

Given a network of stations, a radio-source catalogue, and an interval 
of time of antenna availability, optimize an observation schedule for 
the estimation of, for example, baseline and earth orientation para- 
meters. This problem can al»< 'e approached in two ways. It can be 
asked whether required earth or3 intation variation accuracies mentioned 
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above can be resolved In an olght-hour dally shift rr whether continuous 
observations are needed » the answer being of obvious economic signifi- 
cance. In a similar vein, Mollnder [1978] reported on a method to com- 
pute required antenna time to achieve a given baseline accuracy. 
Alternatively, Ma [1978] searched for a scheduling strategy to minimize 
the baseline component variances. The second-order design problem has 
been somewhat constrained In the past by the low number of sources 
acceptable for geodetic applications (up to 20 with the Mark I system). 
The new Mark III system will enable more sources to be observed 
resulting In a better sky distribution. 

Hints to a possible solution of the design problem are pre" 

sented below. The partial derivatives of the observable with respect 

to a particular parameter constitute the elements of the design matrix 

T 

A that forms the normal matrix N ■ A PA, whose Inverse yields the 
a priori covariance matrix of parameters. In VLBI, the partlals are 
diurnal sinusoids, the baseline vector and source declination deter- 
mining the amplitude, the source right ascension and epoch of observa- 
tion, the phase with respect to o\fT (for example). For second-order 
design, the station locations are given which leaves variable the 
choice of sources and their epochs of observation. The magnitude 
of a particular partial which reflects the sensitivity of an observation 
to a particular parameter determines its numerical contribution to the 
normal matrix. Assumlr lo c;>rreiations between parameters (and obser- 
vations) a diagonal normal matrix would result, and in this case the 
larger the diagonal elements, the smaller the parameter variances. In 
this ideal case, the solution to the design problem would be to observe 


66 


the sources when they would naxlmlxe the partlals with respect to the 
different parameters. However, the parameters are correlated, and Ma 
[1978] found that this approach does not yield optimal results for base- 
line length recovery. Therefore, a particular source must be observed 
not only at the epochs at which the partial derivative sinusoids attain 
their maxima. How then should the observations be distributed? Since 
the partlals enter the normal equations squared, the sinusoids are com- 
posed of two equivalent 12-hour half-cycles. It Is reasonable to assume 
(although correlatlo«i: it«i>tween parameters may Invalidate this assump- 
tion) that observe ^ a both half-cycles are unnecessary from the 
point of view of added sensitivity (though they do add redundancy— on 
the other hand from a systematic error modelling viewpoint, Shapiro 
[1978] suggests observing high declination sources ("clock stars") over 
a large fraction of the diurnal cycle to correct for the effects of 
long-term drifts In the clock behavior) . Observational constraints 
(described below) may limit the availability of a particular source to 
the quarter-cycle. It will be tested by simulations whether the sensi- 
tivity of the observable to a particular parameter can be adequately 
exploited by observing sources throughout a half (or quarter cycle) of 
the corresponding parameter partial. In such a way that the sinusoid Is 
adequately represented. In this way, the sensitivity of the observable 
to a particular parameter may be fully ei^plolted. [Ma, 1978] found that 
the strategy of maximum sources is not optimal. According to the above 
hypothesis this Is due to Inadequate sampling of the sinusoids since 
less observations are available to a particular source while at the same 
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tine antenna slew time Is Increased. Therefore, it would be advan-* 
tageous to observe less sources; Ma suggests ten for geodetic pur- 
poses . 

Until further results are available, the following two-stage 
procedure Is suggested. In the first stage, the sources to be observed 
are selected from the available source list as follows. Suppose we are 
Interested in the optimal estimation of earth orientation parameters. By 
an examination of the partial derivatives It can be seen that for a parti- 
cular station configuration, estimation of the elements of the parameter 
set are sensitive to either low, medium or high declination sources [Bock, 
1980]. Further suppose that It has been decided to observe twelve sources 
over a 24-hour period. Thus, for each parameter we can choose four sources 
whose right ascensions are distributed fairly evenly over 24 hours . These 
sources can be chosen by sorting through the available source list and 
choosing for each group those sources that provide the largest partial deri- 
vative values (a function of source declination) . Once the sources are 
selected, the second stage will Involve choosing the corresponding epochs of 
observation according to the hypothesis v?uggested in the previous para- 
graph . 

There are problems with the above procedure. In a multl-basellne 
experiment there are several baselines to consider. The source sort Is 
then performed according to the "best" baseline for the estimation of a 
particular parameter determined by comparing the sensitivities of the cor- 
responding partial derivatives. In addition, low declination sources are 
visible for shorter periods of time which will require their more Judi- 
cious selection. Finally , and most important, this procedure does not 
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consider the correlations between parameters. Clearly, more research Is 
required into this problem area to translate the above suggestions Into 
mathematical form, or to search for other more rigorous optimization 
techniques. Observation scheduling is quite a tiresome chore and an effi- 
cient algorithm that could provide an ”optlmal" schedule Is needed. 

In developing an optimal schedule algorithm several constraints 
must be considered. A source must be observable from all participating 
stations at a particular epoch of observation. VIP contains a routine 
that outputs visibility matrices for each source over a 2A-hour period 
as described In Section 3.4. An example Is given In Appendix B. 

The slew rates of the station antennas are other factors to be con- 
sidered. Slew rate Is a function of the size, steering mechanism and mount 
geometry of a particular antenna. With equatorial mounts, slew time 
between any two sources Is constant over the entire day, although high 
declination sources are difficult to track. The equatorlally mounted 
Ft. Davis antenna has an hour angle constraint of 5-1/2 hours on each 
side of the meridian. For az-el mounts the slew time between any two 
sources Is a function of the epoch of observation. In addition, there 
Is a blind spot at the zenith as well as cable wrap problems. 

Robertson [private communication] has suggested a function that would 
weigh cost, corresponding to slew time, against benefit to the objec- 
tive function of some optimization technique. 

Another constraint is dead time which is the time required for 
nonobservatlonal matters such as the switching of tapes and water vapor 
radlometry. Economic constraints may Include the availability of only 
one eight-hour shift per day. 
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4.3 Obggrvatlon Correlations 


In Section 3.2.6 a model suitable for covariance analyaea was 
developed for determining the correlations between simultaneously 
observed time delays (and time delay rates) in an N station network. 

The model is a formulation of the theoretical time delay definltion-- 
the difference in times of arrival of a given segment of a wavefront at 
two antennas. However, since delays are not measured in this manner, 
rather by the cross-correlation procedure described in Section 2.3, 
this model needs to be modified for application to real data. Neglect of 
the real observation correlations (or a suitable approximation to them) 
in the second adjustment mentioned in Section 4.1 may alter significantly 
the geodetic parameter estimates. 

In this context, an experiment was performed at the Goddard 
Space Flight Center with the aid of Jim Ryan and Chopo Ma. Two good 
data seta were edited to retain all good simultaneous observations from 
the Haystack, OVRO (Owens Valley) and NRAO (Greenbank) stations. Least 
squares estimates of baseline and earth orientation parameters using all 
three baseline observation sets (a dependent set, see Section 3.2.6) were 
compared to the results of each of the three two-baseline independent 
combinations. Theoretically, the two-baseline combinations should yield 
identical estimates when the true observation correlations are con- 
sidered [Shapiro, private communication]. The three-baseline case 
includes only two independent baselines and therefore will yield overly 
optimistic estimates. In practice, due to inadequate knowledge of 
observation correlations, a diagonal variance- covariance rnatrls: of 
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ob8«rvables 1« assunad. Since the NASA aofcvare cannot acconnodate 
off-diagonal elementa for the variance-covariance matrix of obaerv- 
ablea, the experlmenta were performed with a diagonal matrix. 

The data aeta were divided Into three Intervale of time. The 
flrat aet makes up the flrat Interval of obaervatlona of approximately 
25-hour duration. Thla entire Interval la used for the flrat earth 
orientation atep, thereby providing reference orientation for the polar 
motion and earth rotation parameters of the subsequent steps # as 
described In Section 3.2.3. The second data set la divided Into two 
approximately 20-hour steps. For each of these two steps » three earth 
orientation variation parameters are estimated relative to the Initial 
orientation provided by the first step. 

The experiments were run In the "unweighted" and "weighted" 
observation modes. The unweighted mode Involves the original time delay 
observations with their estimated standard deviations as they are 
recovered from cross-correlation of the tapes. The weighted mode uses 
observations that have been scaled after an Initial adjustment to reduce 
the a posteriori variance of unit weight to unity. Each baseline Is 
scaled differently, the scale factors computed by a numerical procedure 
[Robertson, 1975]. It is felt by those Involved that this scaling 
tends to compensate for unmodelled systematic effects as well as for the 
neglected observation correlations. This writer is not aware of any 
statistical justification for this scaling. 

The results are presented In Tables 4.1 and 4.2. It can be seen 
in both tables that the baseline "components" T, e and a may differ by 
as much as 50 cm from one solution to the next, and In some cases these 
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discrepancies do not fall within the estimated noise levels. Baseline 
distance estimates, on the \)ther hand, are well behaved differing by not 
more than 6.6 cm and are within the noise levels. The polar motion 
variation parameter estimates, AE|£ and are very erratic, while 
AKjj^ variations are well behaved. Of course, in the three baseline 
configuration, the standard deviation estimates are lower than in the 
two-baseline combination cases, but not significantly. However, these 
differences should become more pronounced as the number of stations is 
Increased. Finally, there seems to be little difference compared to 
the discrepancies of the weighted mode, some Increase and some 
decrease. Therefore, at least in this experiment, the "weighting" pro- 
cedure is ineffective in generally redi'cing the discrepancies. 

It is planned to study the effect of Including observation cor- 
relations using these data sets . It is hoped to be able to reduce the 
discrepancies in this manner. However, systeiiiatlc errors may be a more 
important factor in causing these differences than the random nature of 
the observations. In order to test whether correlation neglect is sig- 
nificant, the simplified model of Section 3.2.6 must be at least 
expanded to handle time delays of varying precision. A substantial 
reduction in the discrepancies will Indicate a good correlation model, 
while remaining differences will be due to unmcdelled systematic 
effects. The latter need to be reduced by better instrument calibra- 
tion and by Improved mathematical models. 
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APPENDIX A 


VLB I COVARIANCE AN/XYSIS INTERACTIVE PROGRAM (VIP) 

JCL, Explanatory Tables and Figures, Documented Listing 

A.l Introduction 

In this appendix a documented listing of VIP Is presented, as 
well as the JCL and explanatory tables and figures for the user's ease 
of adaptation and operation. The program, written In FORTRAN, must be 
loaded with the FORTRAN Library (FORTUB), the IBM FORTRAN Scientific 
Subroutine Package (FORTSSP) and the Tektronix Graphics 2 package 
(TXGRAPH2) to achieve Its full capability. The FORTSSP Is called for 
normal matrix inversion (DSINV) and calculation of the normal matrix 
eigenvalues (DEIGEN) . Thus, any other routine that performs the same 
functions may be substituted, though It must be able to handle matrices 
whose upper triangular elements are stored In vector format. The 
graphics portion of VIP may be skipped so that TXGRAPH2 Is optional. 
Consequently, the program may be run on any Time Sharing Option (TSO) 
compatible interactive terminal. 

VIP is mainly intended as a covariance analysis program as 
explained in Chapter 3. However, it is also possible to perform a 
standard least squares estimation of the parameters and their standard 
deviations (a posteriori) but only for simulation purposes. An example 
is given in Appendix B, The program is not equipped to handle real data. 


2 


All the VIF JCL listed In Figure A.l Is given In the form of a 
command procedure (CLIST) of the IBM 0S/VS2 TSO Command Language [IBM, 
1978]. In this case, the program Is stored In a sequential data set 
called BOCK. FORT. The CLIST, stored In BOCKLIB. CLIST, allocates the 
necessary files, compiles BOCK. FORT and loads BOGK.OBJ with 
SYS2.TXGRAPH2, SYSl.FORTLIB and SYS2.F0RTSSP described earlier. The 
contents of each file are listed In Table. A.l. The £..alre procedure Is 
Initiated by the following sequence of commands: 
enter: EXEC BOCKLIB. CLIST 

terminal response: ENTER POSITIONAL PARAMETERS DSNAME 

enter: BOCK. FORT 

It may taken from a few seconds up to a few minutes until the program 
Is compiled and loaded depending on the system status. The first pro- 
gram prompt to the user will be 

DO YOU WISH TO DRAW MAP? 

At the end of the session, the object module, BOCK. OBJ Is deleted. The 
program (2770 card records) and the other necessary data files occupy 
approximately 35 tracks of disk space. 
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Figure A.l VIP CLIST 


Table A.l. VIP File Allocation 


File No. Content 

3 Station file followed by Radio Source fJle 


Format 


(unformatted) 


12 
3A4 

13.12, F6.3 

14.12, F6.3 
F10.2 

12 
3A4 

13.12, F6.3 

14.12, F6.3 

4 Digitized map coordinates 
(optional) 

Format: standard digitizer card format 

6 points per card (8X, 12F6.3) 

5 TSO terminal Input file 

6 TSO terminal output file 

7 Line printer output file (may be VERSATEC) 

8 Card punch file (not used In program) 


first record - 

ellipsoid equatorial radius : . 
Inverse flattening 
for each station record - 
station number 
station name 
latitude (D.M.S.) 
longitude (east) (D.M.S.) 
ellipsoidal height (m) 
for each source record - 
source number 
source name 

right ascension (H.M.S.) 
declination (D.M.S.) 


9 Planned observation schedule (filled prior to run) 

(optional) 

Format: Source number - one per record unformatted - 

to be used only for simultaneous observations 
from all participating stations at even 
Intervals of time 


or Format; Source number, hour and minute of observation - 
one set record unformatted - to be used only 
for simultaneous observations from all 
participating stations at uneven intervals 
of time. 


Table A.l (continued) 


10-15 Simulated observation flies - one per baseline - filled In 
inMonai ^ order of baseline selection, e.g., 10 - first baseline, 

11 - second baseline, etc. 

Format 

each record - 

baseline number 13 

source number 15 

Format 

h 

hour of observation relative to 0 UT 15 

of Initial day 

minute of observation 15 

delay (m) (zero If Interested In F20.10 

covariance analysis only) 

delay rate (m/hr) (zero If interested in F20.10 

covariance analysis only or delay 
observations only) 

index - only for nonslmultaneous observe- 12 
tlons (when not all participating 
stations observe at each epoch) 

0 - next baseline observation at same epoch 

1 - next baseline observation at next epoch 

Use these files for the following cases: 

when entering observation prior to program run, and/or 
when entering a schedule of non-slmultaneous observa- 
tion 


A. 3 Explanatory Information 

Table A. 2 provides • for easy reference, an Index of the various sub- 
routines and their respective purpose,, Flgvre A. 2 Illustrates the flow 
of the program and the Interconnecitlons among the subroutines. 

Table A. 3 contains a listing of all data that needs to be 
Input by the operator at the terminal (and, optionally, prior to the 
program run). The operator Is prompted to supply the Information by 
messages on the screen. Some of the Input may not be requested depending 
on the program options as Indicated In the table. In Table A. 4, the VIP 
program options are listed. These are chosen by the operator Inter- 
actively In response to program prompts. 

There are certain program parameters that may need to be modi- 
fied depending on the user’s needs. These are Indicated at the appro- 
priate locations in the program and are summarized here. In the main 
program, the variables NSTAT and NQUAS refer to the number of stations 
and number of sources respectively, stored on flTe 3 (see Table A.l). 

The values specified in the program are 6 and 47 respectively (see VP 
112, In the listing). A greater number will necessitate Increasing the 
dimensions of the appropriate arrays (see VP 520 - VP 980) . Kemember 
that the maximum available storage on TSO is normally 256K (default 
value - 192K) . 

I 

Subroutines MAPDRW and BSLN assume that the coordinates of the 
United States are digitized, at a scale of approximately 1 ; 10,000,000, 
as well as the station locations of Westford, Owens Valley, Goldstone, 

Ft. Davis, Greenbank and Richmond. Any deviation from these assumptions 
will necessitate minor modifications in these routineiji (see the listing) . 


SO 


Of course » the user must supply a set of digitized map coord’*.nates 
of his area of interest and may need to redefine the screen and virtual 
windows (see VP 1740 - VP 1790) to take Into account the map scale. 
However, the graphics portion of the program Is optional so that 
digitized map coordinates are not a necessity. 

In the Interactive mode the user Inputs data at the terminal 
when prompted to do so by the program. All Input is accepted after the 
RETURN key Is hit. If an error Is made before RETURN, simply hit the 
BREAK key and re-enter. If RETURN has been specified, the program will 
usually provide additional chances, immediately or at a later stage, 
until an acceptable response Is made. However, certain erroneous 
responses will cause the program to abnormally terminate. Therefore, 
it is good practice to examine your responses before hitting RETURN and 
to follow directions carefully. 


81 


Table A. 2 VXP Subroutine Index 


Haae 

Main Program 

MAPDRW 

(optional) 

BSLN 

SIDTIM 

GRESID* 

JULIA* 

STAINS 

QUASAR 

SIMULT 

FLAGS 

WEIGHT 

PARTDR 

AMATR 

FAMTR 


* Adapted 


Purpoae 

Admlnlatera the following: 

MAIN 1: Baseline configuration display 

MAIN 2: Mutual visibility outllner 

MAIN 3: Schedule simulator 
MAIN 4: Least squares estimation 

Plots digitised map coordinates, station locations and 
station symbol selection menu. 

Inputs station and baseline selections and displays them 
on map. 

Inputs time Information and outputs C’'T of Initial epoch 
and chosen Interval of olservatlons. 

Calculates GST of Initial epoch. 

Converts Universal Time to Julian date. 

Inputs station Information and computes baseline coordi- 
nate differences and baseline lengths. 

Sources are displayed and selected. Computes mutual 
visibility matrix* (optional) . 

Simulates oLaervatlons for chosen schedules. 

Inputs experiment flags. 

Inputs observation weighting Information and computes 
weight matrix of observables. 

Calculates partial derivatives of observables with respect 
to parameters. 

Fills design matrix (A) with calculated partial deriva- 
tives for delay and combination of delay and delay rate 
observations. 

Fills design matrix (A) with calculated partial deriva- 
tives for delay rate observations. 


from [Dermanls, 1977]. 
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Table A. 2 (Continued) 


Name 

FILL 

SOLVE 

STDLST 


RAD 

DEGMS 

MATPV 

LOG 

FRAME 

UNITS 

RECT 

EQUITR 

CIRCLE 


Purpoee 

T T 

Fills normal matrix (A PA) and U matrix (A PL) sequen- 
tially. 

Computes variance-covariance matrix of parameters 
(a priori and a posteriori) > parameter correlation matrix 
and normal matrix eigenvalues. 

Computes and outputs estimated standard deviations of 
parameters (a priori and a posteriori) and outputs cor- 
rections to approximate parameters. 

AUXILIARY ROUTINES 

Converts angle In degrees, minutes, seconds, to radians. 

?'%rt'urms opposite function of RAD. 

Performs matrix multiplication for matrices stored In 
general or triangular storage. 

IBM SSP routine-matrix storage manipulator. 

PLOTTING ROUTINES 

Frames a screen window. 

Converts centimeters to virtual coordinates. 

Plots a square. 

Plots an equilateral triangle. 

Plots a circle. 
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Tabla A. 3. VIP Input Paramatara 


Variable Name 

Daacrlptlon 


IN 

Number of atatione 

BSLN 

TNB 

Number of baaellnea 


I8T‘ 

Station aelectlon 


NEX 

Experiment number 


Symbol Selection^ 

Operator moves cursor to choose station 
symbol . 


IYEAR,IMO,IDAY 

IHOUR,IMIN,SEC 

Initial epoch of observations (UT) 

SIDTM 

JYEAR,JMO,JDAY 

JHOUR,JMIN,SECJ 

Final epoch of observations (ITT) 


IQUAS‘ 

Chosen source numbers 

QUASAR 

ZNTMAX^ 

Maximum source zenith distance 


IPFIX 

Reference right ascension source number 


wn 

Number of steps for earth orientation 

MAIN 3 

SFNC^ 

Final epochs of earth orientation steps 


PMX^’*^ 

Approximate step values - Irst 
component of polar motion 'C) 

SIMULT 


Approximate step values second 
component of polar motion (n) 



Approximate step values (OTC-UTl) 


DT 

Time interval between observations 


iFiu;* 

Storage files for observations - one 
per baseline 

(or MAIN 3) 
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Table A. 3 (Continued) 


Variable NaM 

Description 

Subroutine 

K* 

Scheduled Observations - source number 
(when DT Is specified) 


K,IHOUR,MIN* 

Scheduled observations - source nund>er, 
hour (relative to Initial epoch), minute 
(when DT not specified) 


FLAGl, FLAG2 
FLAG3, FLAG4 

Program flags (see Table A. 4 for details) 

FLAGS 

SIGl 

Delay standard error (m) 

WEIGHT 

SIG2 

Delay rate standard error (m/hr) 


Pl‘ 

Covariance matrix of observations - 
upper triangular and diagonal elements 
scaled to unity (one N x N block - see 
Section 3.2.6) 



Prior to ProRram Run^ 


K 

Source number - see Table A.l, file 9 for 

details 

or 



K, IHOUR, MIN 

Source number, hour and minute of observation 

NB, IP, IHOUR, 
IMIN. DS*, 
FRNG®, XEND^ 

Baseline number, source number, hour of observation 
(UT) , minute of observation, delay, delay rate, end 
of observaclon Index - see Table A.l, files 10-15 
for details 


*Arrr 

^Optional. 

®May be set to zero for covariance analysis. 
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Variable 
MAP - 0 

MAP - 1 

GG - YES* 
GG “ NO 

IFY - 0 
IFY - 1 

FG - XES* 

FG - XO 
GG - XES 

GG » XO 

IFLAG - 1 
IFLAG - ?, 

DT » 0 
DT » X 

IPASS « 0 
IFASS « ] 




Table A. 4, VIP Program Optlona 
(Specified by Uaer Interactively) 


Option 

Baseline configuration map is not displayed - 
subroutine MAPDRW is skipped - graphics 
display terminal not required - file 4 ic empty 

Digitized map coordinates from file 4 are 
plotted on graphics screen + baseline configura- 
tion - option to terminate program if only map 
is desired 

Change time input from previous run 
Keep same time input as in previous run 

Skip source mutual visibility outliner - choose 
sources directly 

Compute visibility matrix - plotted source by 
source on screen until specified number of 
sources chosen - can generate complete 
visibility matrix by specifying last source on 
list as last chosen source - may terminate 
program at this point if only visibility matrix 
is desired 

Keep same earth orientation step input from 
previous run 

Change earth orientation step input 

Skip subroutine SIMULT - store schedule infor- 
mation (and optional observa*-ions) on files 
10-15 prior to program run (see Table A.l) - one 
baseline per file 

Call subroutine SIMULT - either store schedule 
information on file 9 previous to run or input 
schedule interactively (see Table A.l) 

Enter observation schedule at terminal 
Input observation schedule from file 9 

Time interval between observations is variable 
Time interval between observations is X minutes 

Observations scheduled ever DT minutes/ 
Observations scheduled at uneven Intervals 


Subroutine 

MAINl 

NAIN2 

QUASAR 

MAIN3 

MAIN3 

SIMULT 

SIMULT 

SIMULT 
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Table A. 4 (Continued) 


Variable 


Option 


Subroutine 


IFRNG - 0 Simulate delay observations only SIMULT 

XFKNG " 1 Simulate delay rate observations, too 


" 0 
ISIM - 1 


FG ^ XES^ 
FG - XO 

FLAGl - 1 
FUGl - 2 
FLAGl « 3 

FUG2 - 1 
FLAG2 “ 2 

FLAG2 - 3 


FUG3 - 1 
FLAG3 - 2 

FLAG4 - 1 
FLAG4 » 2 

FG - XES^ 

FG - XO 

GG - XES 
GG - XO 

ICODE - 1^ 
ICODE “ 2 
ICODE - 3 
ICODE - 4 


All observations are performed simultaneously MAIN4 
from all participating stations 

Opposite of ISIM - 0, when mutual visibility 
makes observations from all stations at a parti- 
cular epoch Impossible 

Keep same flag Input as In previous run MAIN4 

Reinitialize program flags 

Delay observations only FLAGS 

Delay + delay rate observations 
Delay rate only 

Multi-baseline configuration 
One baseline - estimate first component 
of polar motion variations (C) 

One baseline - estimate second component 
of polar motion variations (o) 

Covariance analysis only FLAGS 

Complete least squares estimation 

Estimate all parameters 

Delete clock parameters from parameter list 

Keep same observation weight Input as In previous 
run 

Input new observation weight data 

Rerun program with new data input (see ICODE) 

Terminate session 

Change station input (but not source^ 

Change source input data (but not station) 

Change both station and source input data 
Change other Input data (but not source or station) 


FLAGS 

MAIN4 

MAIN4 


FLAGS 

« 


^For piogtam rerun only 
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A. 4 VIP Documented Listing 

Most of the Information presented In the tables and figures are 
also described In the program documentation. The VIP documentation con- 
sists of a heading at the beginning of each subroutine and other comment 
cards Interspersed throughout the program for added detail. Each 
heading Includes the following Information when relevant; 

1. Subroutine function (title) 

2. INPUT parameters - passes to the routine through the parameter 
list or via common blocks 

3. READ parameters - read within the subroutine using file number 5 

4. WRITE parameters - written from within subroutine using file 
number 6 or 7 

5. OUTPUT parameters - output for use In other parts of program 
by the parameter list or via common blocks 

6. OPTIONS - subroutine options 

7. SUBROUTINES - called by routine 
VIP is listed on the following pages. 
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V»P “ ^ VI-BI INTRRAGTIVR PROGfUH 


PURPOSE I ro AID IN niE SINULATION AND DESIGN OF VLDl 
EXl’ElUMEN'l’S. 

DESCRIPTION t VIP IS INTENDED PlilMAIULY AS A MULTl- 
OPTIONAR GOVAIUANGR ANAI.VSIS PHOGRAH Al.THOUGH IT 
IS POSSIBLE TO PERFORM A LEAST SOUARES ADJUSTfflSNT 
ON SIMULATED DATA. THE PAiUMETER SET INCLUDES 
BASRI.INR VRGWR PARAMETERS, VARIATIONS IN POLAR MOTION 
AND EARTH ROTATION AVERAGED OVER EARTH ORIENTATION 
STEi'S, SOURCE PARAME'rEUS AND CLOCK PAIIAMETEUS (SEE 
CHAPTER 0 FOR AN EXPLANATION OF THE ESTIMABLE PARA- 
METER SET). THE ANALYSIS I’AY INCLUDE DELAY AND/OR 
DELAY RATE OBSERVATIONS. 

VIP IS RUN IN THE INITIRACTIVE MODE ON ANY 
TOO GOrn'ATIULE INTERACTIVE TERMINAL AI/I'IIOUGU TlUfi 
OPTIONAL GRAPH I GS CAPABILITIES ARE DESIGNED FOR 
TEKTRONIX TERMINALS ( E . G . TEKl'RON I X 4012) . IN THIS MODE 
THE USER IS ADLE IG SIMULATE AN EXPEIUMENT, VIEW TUB 
RESULTS IN REAL TIME AND RERUN THROUGH THE PROGRAM 
WITH THE OPTION OF CHANGING ANY OR ALL OF THE PREVIOUS 
INPUT PAUAtUSTERS.TlIIS PROCESS MAY HE lUOPEATED AS MANY 
TIMES AS DESIRED WITH ONE LOADING OF THE PROGRAM. VIP 
PROVIDES EASE OF OPERATION AS WELL AS SAVINGS IN i’lM 
AND COST lUSLATIVK TO THE BA'I’CH MODE. VIP MUST BE 
LOADED WITH FORTLIB, FOHTSSP AND TXGRAPH2. 

VIP IS DIVIDED INTO FOUR SECTIONS ! 

t. BASELINE CONFIGURATION DISPLAY 

2. MUTUAL VISIBILITY OUTLINEU 

3. SCHEDULE SIMULATOR 

4 . LEAST 80UAHE8 ESTIMATION 

WRITTEN BY YEHUDA DOCK DEPT GEODETIC SCIENCE 1979 
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** 

tut 
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MAIN PROCUAJI 


IMPLICIT REAL*R( A-H,«-Z> 

REALtf'4 XX.YY.ARAY 
lN'l'EGERiH2 INDEX 

INTEGER TNB.FI.AGI ,FLAG2,FLAG,T,FLAG4 

MAXIMUM NUMUEU OF PAIiuAMETERS ALLOWED HY DIMENSIONS 

1. ft BASEUNES 

2. 12 OUASARS 

U, 4 KAK’I’R ORIENTATION S'l'EPS 
4. ft CLOCK RATE PARAMETERS 

3 . 6 CLOCK OFFSET PARAMETERS 

TIHC 'WTAL NUNRIOl OF PAHAME'I'EHS SHOULD NOT EXCEED ft2 
INCKEASE DIMENSIONS FOR A LARGER PARAME’I’ER SET 
NOTEiflAXlMUM STORAGE AVAILADLE ON TOO-206K 


ARRAY FUNCTION 

ARAV TERMINAl. STATUS ARRAY 

XX, YY DIGITIZED MAP GOORDINA'IES 

1S,JS,1ST INDICES FOR CHOSEN STATIONS 


SUBROUTINES 

MAPDRW.BSLN 
MAPDRW, DSLN 
RSI.N.STATNS 
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90 
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110 
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I HO 
140 
150 
tftO 
170 

tuo 

190 
200 
210 
220 
200 
240 
250 
200 
270 
280 
290 
HOO 
310 
020 
330 
340 
000 
3<i0 
070 
OBO 
390 
400 
410 
420 
400 
440 
430 
4ft0 
470 
4U0 
490 
500 
510 
520 
530 
040 
550 
560 
570 
5B0 
690 
ftOO 
ft 10 
ft20 
ftno 
040 
630 
ftftO 
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G 

C 

G 

0 

G 

C 

0 

c 

c 

G 

C 

C 

G 

C 

G 

G 

G 

G 

G 

C 

0 

G 

G 

G 

G 


G 


,7.T 


X.Y.Z 
Kl’, YV, 

DIST 
XG, VG.ZG 
XI’I>,VTP,7.TP 
im, DEGR 
El.EU.EU 
KT.KP 
iKDEX 
IttUAK 
RA.D 
IFILE 

l'MX,HMY,PW7 

AL 

IM'l.PU 
PART.G 
AM, AG 
UG 
W 

u 

XA 

B.DM.EM 

GUUU 

8KHG 


CAUTEBIAN STATION GOOUUINATES 
STATION UNIT VEGWR GOMI'ONENTS 
BASELINE DISTANCES 
BASELINE GOOIUUNATE UIPPEIIENGES 
GHOSEN STATION UNIT VECTOR CONPON. 
QUASAR COORDINATES ON PILE 
QUASAR UNI'l' VECWm GOMPONKN'I'S 
TITLE FOR OBSERVABILITY OU1TINER 
OBSERVABILITY MATRIX 
INDICES FOR CHOSEN QUASARS 
CHOSEN QUASAR COORDINATES 
OUTPUT FILES FOR OBSERVATIONS 
APPROX VALUES EARTH ORIENTATION 
MISCLOSURE VECTOR 
WEIGHT WA'l’RIX MAN IPULA'IORS 
PARTIAL DERIVATIVES DELAY ft RATE 
DESIGN MATRIX MAN I PULAIORS 
WORK MATRIX (BC=ATP) 

NORMAL MATRIX AND ITS INVERSE 
U=ATPA VECTOR 

CORRECTIONS I'D APPROX PARAMETERS 
PROPAGATION OF BASELINE DISTANCE 
GOllRELATION MATRIX MANIPULATION 
STEP FINAL EPOCHS 


STATHS 

STATNS 

STATNS, SOLVE 

STATNS , S I MULT, PART 

STATNS, QUASAR 

QUASAR 

QUASAR 

QUASAR 

QUASAR 

QUASAR 

QUASAR, S IMULT, PART 

SI MULT 

SI MULT 

PARTDR,FILL 

WEICHT,FILL 

PARTDR, AMATR 

AMATR,FILL 

FILL 

FILL, SOLVE 
FILL, SOLVE 
SOLVE 
SOLVE 
SOLVE 
MAI NO 


DIMENSION RA(12> , D( 12) , KEEP(24), IQUASI 12) , INDEX(6,24), AIUY(60 
I). Xdi), Vlft), 7,(R). IS(<)). .18(1)) . ISTdi), IFILEUO. RAIU47) . DECR 
2(47), EK47), E2(4T), £0(47), SIHO), YT(w), 7T(6) , PMX(4), PMY(4>, 
DPMZ(4), XC(6), YC((>), ZC(0) , DIST(ft), W( I45!l) , U(ft2), XA((i2), PAR 
4T(I0). G(l«). AM(124). G0RR(<.2), B( 100) , BM(0(i), EM( 108) , XTP(«) , 
5VTP(6), ZTP(6), AG(744), DG(744) , !*,7B) , PU(2(), XX(6), VY(6), Pl( 
642), AL( 12> , SFNC(4) 

COMMON /DRAWl/ XX/DRAW2/YY/DRAWa/ARAY/CRDl/X'l’l' '0RD2/Y'1’P/G1U)U/7TP/C 
I lU)4/XC/CRI)n/YC/CRDf./7.C/Tl ME/ 1 NO , I DAY, I VEAIVDEXI / IS/DEXa/.FS/DEXfl/ IS 
2IVDEX4/NEX/S0U 1 /RA/S0U2/D/FLG/FLAG 1 , FLAG2 , FLAGO , FLAG4/PDR i/PAIlT/PD 
BR2/G/ ATllX/ Ari/NT’RX/W/U’l’RX/U/BS/ TNB , KK 


DEFINE TERMINAL RESPONSES 
DATA X0,XES,XEHM,X1LE/*N0* 


'YES’ 


TERM’ , ’FILE’/ 


DEFINE PROCllAM CONSl’ANTS ! 

NS'l’A'P “ TOTAL NUMRER OF STA'ITONS ON FILE (GRANGE IF NECESSARY) 
NQUAS - NUMHF.rt OF QUASARS ON FILE (CHANGE IF NECESSARY) 

IF INCREASING NSTAT AND/OR NQUAS CHANGE Al’PROl’RlATE DIMENSIONS 
I) A TA NSTA'IV (. / . NQU AK/47/ 

OONV “ CONVERSION FACTOR UNIVERSAL TO SIDEREAL TIME 
ERAD “ approximate EARTH RADIUS 
C - THE SPEED OF LIGHT IN NETERS/NSF.G 

DATA CONV/1 .OO2737909265D0/,ERAD/<»a7lDOO.DO/,C/O.29O7D«/ 
l’I»4.D0.tDATAN( l.UO) 

no i CONVERSION FAGIOR FROM RADIANS TO SECONDS OP ARC 
CALL RAD (0,0, 1.0,P,P1) 

U0= I.UO^USIN(P) 

OMEGA IS rilE ROTATION RATE OF THE EARTH IN HADIANS/HOUR 
OMG»PI!liUONV/12.D0 


INITIALIZE 
CALL INIIT 


TERMINAL 
( 120 ) 


CONTROL SYSTEM 


MAINl ! BASELINE CONFIGURATION DISPLAY 

GCGCQGCCCCGCGCCGGCCCCCCCOCCCUGCCOCCCCCGCGCGGCGGCGCGCCCGCGGCCCGGUC 

PURPOSE t CHOOSE PARTICIPATING STATIONS AND PROVIDE 
GRAPHICAL DISPLAY (MAP) OF EACH EXPERIMENT. 

MAP IS OPTIONAL SO THAT THIS SECTION MAY 
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670 
680 
690 
700 
710 
720 
730 
740 
700 
760 
770 
7110 
790 
800 
BIO 
820 
UU0 
840 
800 
860 
870 
800 
1190 
900 
910 
920 
mii 
940 
•>50 
960 
970 
980 
990 
1000 
1010 
1020 
1030 
1040 
1000 
1060 
1070 
1000 
1090 
1100 
1110 
I 120 
1130 
1 140 
1100 
1160 
1 170 
1180 
1 190 
1200 
1210 
1220 
1230 
1240 
1200 
1360 
1270 
12110 
1290 
1300 
1310 
1320 
I ‘330 
1340 
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0 
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OT, RUN ON /iNV INTF.RACTIVR SCRRRN OR TRRNINAL 
AT mCG AND NOT NECESSARILY A TERMINAL WITH 
CHATUICS CAI’AUILITY. 

OI>riON 10 TERMINATE PROORAM AT END OF SECTION 
WHEN MAP DISPLAY IS ONLY DESIRED OUJECTIVE. 

SEE INDIVIDUAL SUBROUTINES FOR MORE INFO. 

GALLS SUUROUTINES MAPDRW.DSLN 

CCCCCGGCCCCGCCCCCCCCCCCCCCCGCCCCCCCCCCCCCCCCCCCCCCCCCGCCCCCCCCCCC 

IG0DE=0 DENOTES INITIAL PROGRAM RUN (OTHER I CODES ARE PROGRAM 
RERUN OPTIONS - EW*LAINED AT END OF MAIN PROGIUM) 

IC0DE=0 


OPTION ! MAP“0,MAP IS NOT DRAWN 
flAl'=l,MAl' IS DIUWN 

MAP»0 

ERASE SCIU;EN 
GAU, ERASE 

MOVE CURSOR TO UPPER LEFT-HAND CORNER OF SCREEN 
CALL HOME 

SAVE CURRENT STATUS OF TERMINAL CONTROL AREA 
CALL SVSTAT ( ARAY) 

SWnCil lO ALPHANUMERIC NODE 
CALL ANMODE 
WRITE (0,3) 

I'ORMAT ( ’ DO YOU WISH lO DRAW MAI") 

READ (5,4) F 

IF(I' .NE.XES.AND.F.NE.XO) GO TO 2 
EORMA’I' (A4) 

CALL ERASE 
CALL HOME 

lEd’.EO. xo) no 10 .5 

FILE 4 CONTAINS DIGITIZED MAP COORD I NATES- MAY HE EMPTY IF MAP»0 
1F( ICODE.NE.O) REWIND 4 
MAP= 1 

CALL RES'l’AT (ARAY) 

DEFINE SCREEN AND VIRTUAL WINDOWS 
GALL SWINDO (350,000,270,605) 

VIRTUAL NAPS TO SCREEN IN A 2s 1 RATIO - DIGITIZED MAP TO SCREEN 
SCAT.E CAN DE CHANCED TO SUIT PARTICULAR MAP 80 IT FALLS IN 
SCREEN WINDOW UY SUITABLE CHOICE OF VUITUAL WINDOW 
CALL VWINIH) (0.0, 1320. ,0.0, 1010. ) 

CALL MAI* llANDLEll 
CALI, MAPDRW ( HSTAT) 

CALL BASELINE lANDLEU 

CALL nSLN ( 1N,MAP,NSTAT, ICODE) 

CALL ERASE 
CALL HOME 
CALL ANMODE 
WHITE (0,7) 

FORMAT ( • DO YOU WISH TO RUN OR RE-RUN MAP DRAWINC^ SESSION’/’ 
lOU STATION AND UASEHNE SELECTION’) 

HEAD (5,4) no 

IF(GG.NE.Xir.S.AND.GC.NE.XO) GO TO 0 

IF(CC.Ett.XO) CO TO U 

CALL ERASE 

REWIND 4 

CO 10 1 

WRITE (6,0) 

FORMAT (’ DO YOU WISH TO TERMINATE SESSION’) 

READ (5,4) GG 


VP 1850 
VI’ 1300 
VP 1370 
VP 1300 
VP 1300 
VP 1400 
VP 1410 
VP 1420 
VP 1430 
VP 1440 
VP 1450 
VP 1400 
VP 1470 
VP 1400 
VP 1400 
VP 1500 
VP 1510 
VP 1520 
VP 1530 
VP 1540 
VP 1550 
VP 1500 
VP 1570 
VP 1500 
VP 1500 
VP 1000 
VP 1610 
VP 1020 
VP 1030 
VP 1640 
VP 1050 
VI’ 1600 
VP 1070 
VP 1000 
VP 1600 
VP 1700 
VP 1710 
VP 1720 
VP 1730 
VP 1740 
VP 1750 
VP 1700 
VP 1770 
VP 1700 
VP 1700 
VP 1000 
VP lino 
VP 1020 
VP 1030 
VP 1040 
VP 1050 
VP 1060 
VP 1070 
VP 1000 
VP 1000 
VP 1000 
VP 1010 
VI* 1020 
VP 1030 
VP 1040 
VP 1050 
VP 1060 
VP 1070 
VP 1000 
VP 1000 
VP 2000 
VP 2010 
VP 2020 



PAnr. 
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11 
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ir((!O.Nr,.XER.AKD.Gfl.KF..XO) CO I'D B 
IF(GG.EO.XES) GO TO 74 

HAIN2 ! MUTUAL VISIBILITY OUTLINER 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCGCCCCCCCCGCCCCCCCCCCCCCCCCCCCCC 

PURPOSE ! INPUT PARAMR’I'ERS FOR CALCIR.ATION OF MimiAL 
VISIBILITY OF OUASAR FROM ALL PAHTICIPATIMO 
STATIONS. VISIBILITY MATRIX OPTIONAL BUT INPUT 
PABAMETEBS NECESSARY FOR REMAINDER OF RUN. 

IF VISIBILITY MATRIX THE ONLY OBJECTIVE TIDSRE 
IS OPTION OF TERMINATING PROGRAM. 

SEE INDIVIDUAL SUBROUTINES FOR MORE INFO. 

CALLS SUBROUTINES S I DTM.STATNS, OUASAR 

CCCCCCCGCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

CALL ERASE 
CALL ROME 
CALL ANMODE 

OPTION lO RETAIN SAME TIME INPUT FOR PROGRAM RE-RUN 
IF< ICODE.EO.O) GO TO 12 
VRITE (6,11) 

FORMAT (’ DO YOU WISH lO CHANGE INTERVAL OF OBSERVATIONS') 

READ (5.4) GG 

IF(GG.NE.KES.AND.GG.NE.XO) GO lO 10 
lF(CC.EO. XO) GO TO 13 
CALL TIME HANDLER 
CALL SIDTM (TilO.Pl ,TF) 

CALL STATION HANDLER 

CALL STATNS ( X, Y,Z, XT, Yl'.ZT, DIST, PI , ICODE, N8TAT) 

OPTION TO SKIP OUASAR SELECTION FOR PROGRAM RE-RUN 
1F( ICODK.EO. 1) GO IX) 16 
CALL ERASE 
CALL ROME 

CALL OUASAR HANDLER 

CALL OUASAR ( XO, UECR, BAR, El , E2, EU, OMC, PI , IttUAS, KIOEP, INDEX, ICODE, IM 
1, IPFIX.NOUAS.TRO.TF, IN) 

OPTION *1X) lERMlNATE SESSION IF ONLY INTERESTED 
IN VISIBILITY OUTLINER 
WRITE (6,9) 
lUOAD (5.4) CG 

ir(GG.NE.XES.ANI).GG.NE.XO) GO TO 15 
IF(GG.EO.XES) GO IX) 76 

MAIN3 s SCHEDULE SIMULATION 

CCCCCCCGCGCCCCCCCCCCCCCCCGCCCCCCCCCCCCCCCCCCCCCCCCCGCGCCCGCCCCCCC 

PURPOSE ! SimiLATE CHOSEN VLB I OBSERVATION S(.HEDULES 

TIME DELAY AND/OR TIME DELAY RATES MAY BE SIMULATED 
VARIOUS OPTIONS AVAILABLE- SEE SUBROmTNE 
SIMULT FOR DETAILS. OPTION IX) SKIP OBSERVATION 
SIMULATION IF WORKING WITH SIMULATED DATA - IN THAT 
CASE OBSERVATION FILES FILLED PRIOR IX) RUN. 

INPUT t ISTT STEPS EARTH ORIENTATION STEP FUNCTION 
SFNC FINAL EPOCHS OF EARTH ORIENTATION STEPS 


VP 

VI' 

VP 

VP 

VP 

VP 

VI' 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VI' 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 


2030 
2040 
2050 
2060 
2070 
2000 
2090 
2100 
21 to 
2120 
2130 
2140 
2150 
2HiO 
2170 
21 BO 
2190 
2200 
2210 
2220 
2230 
2240 
2250 
2260 
2270 
2200 
2290 
2300 
2310 
2320 
2330 
2340 
2350 
2360 
2370 
2300 
2390 
2400 
2410 
2420 
2430 
2440 
2460 
2460 
2470 
2400 
2490 
2500 
2310 
2520 
250© 
2540 
2550 
2560 
2570 
2500 
2590 
2600 
26 10 
2620 
2630 
2640 
2650 
266© 
2670 
2600 
2690 
2700 
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c 

c 

c 

c 

c 

c 


c 

c 


17 

lU 


19 


C 

c 


ae 

c 

21 


u 


22 


22 

24 

28 


26 

C 

C 

27 

28 


C 

C 


n 

29 

28 


21 


C 


WUTR ! ISIT.SFWC 

CALLS SUUHOUTIWE 81WULT 


CCCCCCCCCCCCCCCCCCCCCCGCCCCCCGCCCCCCCCCCCCCCCCCCCCCCGCCCCCCCCCCCC 


GALL RRASR 
CALL HOME 
CALL ANNODE 


OE'nOW IX) SKIP STEP INPUT - l<’OU PllOGlUH UE-HUN 
IF( ICODR.RO.O) CO TO 20 
MUTE (6,18) 

FOHMAT (’ UO VOU WISH TO SKIP STEP INPUT’) 

READ (8.4) FC 

IF(FG.NE.XES,AND.FO.NE.XC GO IX) |7 
IF(FC.Ea.XO) UO TO 20 
MUTE (7.19) 

FOlWlAT (1111) 

CO TO 27 


INPUT -»^ST’EP8 AND FINAL EPOCllS FOll EAllTU OIUENTATION STEl’ FUNCTION 
WHITE (6.21) 

CHANCE FORMAT IF MAXIMUM STEP #(4) HAS BEEN INCREASED 
FORMAT (/,’ ENTER #STEPS FOll EARTU ORIENTATION’,/,’ MAXIMUM NUMBER 
1 OF .steps is 4*) 

READ (5.:m) isrr 

CHANCE IF MAXIMUM STEP *(4) UAS BEEN INCREASED 
IF( LSTT.CT.4) GO TO 20 


WRITE (6.22) 

FORMA'I' (/,’ ENTER END OF STEP EPOCH FOR EACH .STEP’/.’ FORMAT ! HOU 
IR, MINUTE ( INTEGERS) •/, • PRESS RETURN THEN ENTER NEXT FINAL EPOCH 
2’/,’ NOTE! HOUR RELATIVE IX) INITIAL EPOCH’/,’ (MAY BE CREATE 

2R THAN 24) ' ) 

WRITE (7.28) ISTT 

FORMAT ( lUl, UK, ’NUMBER STEPS FOR EARTU ORIENTATION PARAMETERS »’,I 
12 ) 

WRITE (7,24) 

FOHMAT (//,UX, ’ STEP FINAL EPOCHS - RELATIVE TX) INITIAL EPOCH’) 

DO 26 1=1. LS1T 

READ (6,Xt) lUOUR.IMIN 

WRITE (7,28) I, UlOUIl, IMIN 

F0RJ1AT (/.10X, 'STEP 12, 4X, 13. 2X, 'HOURS’ ,8X, 13, 2X, ’MINUTES’ ) 

SFNC( I)=DFLOAT( lUOUlU + lMlN/60. DO 
CONTINUE 


OPTION TX) SKIP SCHEDULE SIMULATION 
WRITE (6.2B) 

FORMAT ( ’ DO YOU WISH TO SKIP SCHEDULE SIMULATION ? ’/’ ANSWEIl Y 

lES OH NO’) 

READ (8,4) CO 

1F(GG.NE.X1;S.AND.GG.NE.X0) GO TX> 27 
1F(CC.EU..XES) CO TX) 29 


CALL SCHEDULE S1MULATX)11 

CAI.L SIMULT (Xl’.RM.XILR.XES.ERAD.PMX.PMY.PMZ.OMO.TF.THO, IFILE, ISTT, 
ISFNC. ICODE.CONV) 

GO TX) 22 


WRITE (6,30) 

FORMAT (’ ENTER OUSERVATION FILE NUMUKRS’/’ AVAILABLE FILE NUMBERS 
1 10-18’) 

READ (5, NO ( IFILE(J) ,J» l.’l’HB) 

DO ’Jl .l=l.TNH 

1F( IFILE(.I) .LT. 10.0H, IFILE(J) .OT. 18) GO TO 29 
CONTINUE 

IF( ICODE.F.U.O) CO 1X) 24 


VP 27 !• 
Vr 2720 
VP 2720 
VP 2740 
VP 2780 
VP 2760 
VP 2770 
VP 27U0 
VP 2790 
VP 2200 
VP 21)18 
VP 2820 
VP 2030 
VP 2840 
VP 2050 
VP 2U60 
VP 2870 
VP 2000 
VP 2090 
VF 2900 
VP 2910 
VP 2920 
VP 2930 
VP 2940 
VP 2980 
VP 2900 
VP 2970 
VP 2900 
VP 2990 
VI' 3000 
VP 3010 
VP 3020 
VP 3030 
VP 3040 
VP 0050 
VP 3060 
VP 3070 
VP 3000 
VP 3090 
VP 3100 
VP 3110 
VP 3120 
VP 3130 
VP 3140 
VP 3150 
VP 3160 
VP 3170 
VP 3100 
VP 3190 
VP 3200 
VP 3210 
VP 3220 
VP 3230 
VP 3240 
VP 3250 
VP 3260 
VP 3270 
VP 32110 
VP 3290 
VP 3300 
VP 3310 
VP 3320 
VP 3330 
VP 3340 
VP 3330 
VP 3360 
VP 3870 
VP 33U0 
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f5« na nd on oooft cftnaonoacicineftoftoo 


C 

38 

83 


40 

41 

42 


44 

C 

C 


RRWIfID ORRF.nVATION FILES IF HE-RUWNINC PROGRAM V? 

DO aa I»1,TNB VP 

NUM- IFM.r.l 1) VP 

REWIND NUM VP 

CONTINUE VP 

VP 

MAIN4 > LEAST SOUARES ESTIMATION HANDLER VP 

VP 

VP 


CCCCCCCCCOCCGCCGCCCCCCCCCCCGCCCCCCCCCCGCCCCCCCCCCCCCCCCCCCCCCCCCC VP 


VP 

PURPOSE ! A) COMPUTES NUMBER OF PARAMETERS FOR ADJUSTMENT VP 

D) ZEROS OUT WORK VECI'OHS VP 

C) INPU1S OBSERVATIONS FROM STOIUCE FIIJ5S VP 

0) CALLS LEAST SQUARES ROUTINES VP 

E) PRESENTS PROGRAM RERUN OP’PIONS VP 

VP 

CALLS SUBROUTINES FLAGS, WEIGHT, PARTDR.FAMTH, VP 

AMAI'R, FILL, SOLVE VP 

VP 

CCCCCCGCGCCCCCCCCCGGCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCGCCGCCCCCC VP 

VP 

CALI, ERASE VP 

CALL HOME VP 

GALL ANMOIIF. VP 

VP 

SIMULTANEOUS OBSERVATIONS FROM ALL STATIONS OPTION VP 

iSIM=0 i SIMULTANEOUS OBSERVATIONS VP 

ISIM«1 s NON-SIfRJLTANEOUS VP 

lSm=0 VP 


WRITE (ft.aO) VP 
FOIWIAT (• ARE ALL STATIONS INVOLVED AT EACH EPOCH OF OHSEHVATION' ) VP 

VP 
VP 
VP 
VP 
VP 
VP 
VP 
VP 
VP 
VP 


VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 


READ (5,4) FG 
IFl FG.NK.. WS.AND.FC.NE.XD) GO IV) 35 
IF(FG.EQ.KO) lSm*l 

OPTION TO SKIP FLAG HANDLER » FOR PROGRAM RERUN 
IF( ICODE.EQ.O) GO IX) 3‘> 

WRITE (0,310 

FORMAT ( • DO YOU WISH TO SKIP FLAG HANDLER’) 

READ (5,4) I'Q 

I F( I’C . N F, . JO'S , AND . FC . NE . XO) GO IX) 37 
IF(FG.Ett.XES) GO TO 40 

CALI. PROGRAM FLAGS HANDLER 
CALL FLAGS (JK) 

OPTION TO SKIP WEIGH!’ HANDLER - FOR PROGRAM RERUN 
IF( ICODE.Ett.O) GO IX) 40 
WRITE (0,42) 

FORMAT ( • DO YOU WISH TO SKIP WEIGHT HANDLER* ) 
READ (3,4) FG 

IF(FC.NE.XFS.AND.FC.NE.XO) GO !X) 41 
IF(FG.E(i.XE8) GO TO 44 


CALI. OHSERVATION WEICHT HANDLER 
CALL WEIGHT (P.PD.SIG.Pl) 

COMPin’F, THE NUMBER OF PARAMETERS TO BE ADJUSTED 
KK IS THE NUMBER OF PARAMETERS 
KK=5*TNB+2^ IPI+( 3~.IK)4( IS'IT- 0-1 
IF DELAY HAI'ES ONLY CHANGE NUMBER OF PARAMETERS 
IF( FLAGl . Ett. 3) KK=!jri'NU+2* IM+( 3-JK) *( ISn- 1) -2 
IF NO CLOCK PARAMETERS CHANGE NUMBER OF PARAMETERS 
IF( FLAGl . NE. 3 . AND. FLAG4. EQ. 2) KK=» KK-2S'TKD 
1F(FLAGI .EQ.U. AND.FLAG4.EU.2) KK»KK-TNB 


VP 

VI’ 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 

VP 


KI0i»ICK*2 VP 

K2«KK»(KK+l)/2 VP 


3399 

3499 
3410 
3420 
3430 
3440 
3459 

3400 
3470 
34U0 
3490 

3500 
3010 
3520 
3530 
3540 
3550 
3500 
3570 
3580 
3590 
3000 
3010 
3020 
3030 
3040 
3050 
3000 
3070 
3080 
3090 
3700 
37)0 
3720 
3730 
3740 
3750 
3700 
3770 
3700 
3790 
3800 
3810 
3820 
3030 
3840 
3850 
3800 
3870 
3800 
3090 
3900 
3910 
3920 
3930 
3940 
3950 
3900 
3970 
3980 
3990 
4000 
4010 
4020 
4030 
4040 
4050 
4000 
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«CS»KK 

VP 

4070 


IF(FLAG1.EQ.2) KS«KK2 

VP 

4080 

c 


VP 

4090 

c 

7.F4R0 OUT WORK VECTORS 

VP 

4100 


DO 45 1.K2 

VP 

4110 


W< J)«0.UO 

VP 

4120 

45 

CONTINUE 

VP 

4100 


DO 46 J> l.KK 

VP 

4140 


IH.I)«0.n0 

VP 

4150 


XA( J)«0,D0 

VP 

4160 

46 

CONTINUE 

VP 

4170 


no 47 1« > ,KK2 

VP 

4100 

47 

AMI I)>»O.DO 

VP 

4190 


DO 40 .1= 1 . 10 

VP 

4200 


AI,( .i)»o.no 

VP 

4210 


GliJMO.DO 

VP 

4220 

40 

l'ART(.l)=0.nO 

VP 

4280 


VTS>VI«0.D0 

VI’ 

4240 


lCOUNT-0 

VP 

4250 


IF, NINO 

VP 

4260 


L« 1 

VP 

4270 

U 


VP 

4280 


IF( IGOnE.F.0.4) GO TO 51 

VP 

4290 


WRITE (7,49) 

VP 

4000 

49 

FORMAT < nil, 12X, 'OUSEIIVATION SCURDULEV) 

VP 

4010 


WRITE (7,50) I, 

VP 

4020 

50 

FORMAT ( /,9X, ‘STEl’ ’ , I2/6X, ’USLH’ , IX, ’UUAS’ ,2X, ’lUl’ ,2X, ’MIN’ ,5X, ’ 

VP 

4000 


IDF,?, AY (M) MIX, ’DELAY RATE (M/HR)V> 

VP 

4040 

C 


VP 

4050 

C 

INPyr OUSERVAITONS 

VP 

4060 

C 

1, ! BASF,!, 1 NR COUNTER 

VP 

4370 

c 

IJ s EARTH ORIENTATION STEP NUMBER COUNTER 

VP 

4080 

51 

iJ=0 

VP 

4890 

52 

I.ls 1.1+ 1 

VP 

4400 

C 

UPDATE END OF STEP EPOCH 

VP 

4410 


ST=SFNG( 1.1) 

VP 

4420 


IF( I.I.OT. 1) GO TO 61 

VP 

4430 

C 

READ OUSERVATIONS-DASELINE BY BASEL IN E.EPOGU BY EPUGU 

VP 

4440 

53 

NUM=IFILK(l.) 

VP 

4450 

C 

NB ! BASELINE NUMBER 

VP 

4460 

G 

IP 1 OUASAR NUMUER 

VP 

4470 

G 

IIIOUR, IMIN ! EPOCH OF OBSERVATION 

Vr 

4480 

G 

DS ! DELAY OBSERVABLE (PATH DIFFERENCE) 

VP 

4490 

G 

FRNC ! DELAY RAl'E ODSERVADLE (PATH DIFFERENCE lUTE) 

VI> 

4500 


IF( ISIM.EO. 1) CO TO 55 

VP 

4510 


READ (NUM,54,END=69) ND, IP, IHOUR, IMIN.DS.FRNG 

VP 

4520 

54 

FORMAT (415,2F20. 10) 

VP 

4530 


CO TO 57 

VP 

4540 

G 

lEND ! END OF EPOCH INDEX (NON-SI MULTAN EOUS OBSERVATIONS) 
READ (NUN,56,RNI)=69) NR, IP, IHOUR, IMIN, DS, FRNC, IKND 

VP 

4550 

55 

VP 

4560 

56 

FORMAT (4I5,2F20, 10. 12) 

VP 

4370 

G 

INCREASE OBSERVATION COUN'PEU BY 1 

VP 

4580 

57 

ICOUN’I'^ ICOUNT+ 1 

VP 

4590 


TK» DFLOATI I HOUR) + 1 MI N/60 , DO 

VP 

4600 


1F( lCODE.Ea.4) CO TO 60 

VI’ 

4610 


IFCTK.LR.ST) CO TO 58 

VP 

4620 


UK* IJ. 1 

VP 

4630 


WRITE (7,50) UK 

VP 

4640 

58 

WRITE (7,59) NB, IP, IHOUR, IMIN, DS, FRNC 

VP 

4650 

59 

FOHMT (4X,4I5,2F20. 10) 

VP 

4660 

C 

UPDATE STEP HUMBER IF NECF«SARV 

VP 

4670 

60 

IF(TK.CT.ST) 00 TO 52 

VP 

4680 

G 


VP 

4690 

C 

CAI,CUI,ATE PARTIAL DERIVATIVES FOR PRESENT OBSERVATION 

VP 

4700 

61 

CALL PARTDR ( IP , TK, XC( NB) , YC( ND> , ZG( NB) . ’1110 , OHO, 0 , UO. DS , FRNC.GONV. 

VP 

4710 


lERAD.AL, JK. U,NU) 

VP 

4720 

G 


VP 

4730 

C 

FILL DESIGN MA'llUX WITH PRESENT OBSERVATION PARTIALS 

VP 

4740 
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sr'jg5^®»s^^!3g!SSJ^5!SS8Er- . 


C 

C 


«a 

68 


C 

C 


64 

65 


C 

G 


G 


66 


C 

C 


6T 

68 

G 

G 

G 

69 


G 

C 

G 

G 

C 

G 


70 

71 


72 

73 


74 

75 


G 

C 


76 


fFCFLACl.NE.a) GO TO 62 

FILL DELAV HATE UE8IGN NATIIIX 

GALL FAMTIl ( Nil, IP, IM. M,I8'IT,.IK, IFFIXI 

GO TO 60 

FILL OELAV OH OELAy-OELAV HATE UE8IGN NATIIIX 
GALL AMATH ( NB. IP. IN. LI, ISTT.JK, IPFIX) 

DO 64 IZ>I.KS 

FILL PHEHF.NT EPOGH I'OHTION OF "A" NATIIIX WITH ONE BARELITIE 
GONTIllRimON AT A TIME 
JZ«(ND-I)>KKH+IZ 
AG( JZI-AIK IZ) 

GONTINIIE 
DO 65 KJN»I,KK2 
AM( KN)«O.DO 

CIWGX FUH NEXI' EPOGH OF 0U8EIIVATI0N 

IF( IF.ND.F.O. I) GO TO 66 

L-L+l 

IFIL.LE.TMD) GO 'lO 00 

HE" INITIALIZE BASELINE COUNTER 

L«I 

ADD CONTRIBUTION OF AN EPOCH OF OBSERVATIONS TO HORHAL MATRIX 
CALL FILL (AC,P,VrPVl,AL,UG) 

K*/,*TNH*KK 
DO 67 KN“1,KZ 
AC(KN>»O.DO 
DO 6B KN» I . TNR 
AL(KN)*O.D0 

GON'ITNUE lO WEXr EPOCH OF OBSERVATIONS 
GO TO 5.T 

CAM, HOUJ'I'ION SUimOUTINE 

CALL SOLVE (SIG.CORR.XA, IN, ISrr.DIST.B.DH.EN.VTPVI, ICOUHT) 

CALL ERASE 
CALL HOME 
CALL ANNODE 

PHOGR.IM RERUN OPTION 

1, JCOUE-1 CHANCE HAHELINE GONFIGUIUTION (OR TIME INPUT) 

2. IGODE«=a CHANGE OUASAR SJ;i,ECTION 
a, ICODE»a CHANCE DOTH OF 'HE ABOVE 

4. ICODE-4 CHANGE OIlUiR INPUT PAHAME'I’ERS 
WRITE (6,71) 

FORMAT (/' DO YOU WISH lO HUN niE PHOCHAK AGAIN’) 

READ (5,4) CC 

IF(GG.NE.Xl.S.AND.GC.NE.XO) GO TO 7© 

IF(CC.E(i.KO) CO 'lO 76 
WRITE (6.70) 

FORMT (//• DO YOU WISH PO CHANCE I’llE BASELINE CONFIGURATION’) 

lUl'AD (5,4) CC 

lE(CC.NE.XES.AND.GG.NE.XO) GO 10 72 
iF(cc.Ea.xo) CO ro 74 

1C0DE= I 
WRI'I'E (6.75) 

FOUMAT (//’ DO YOU WISH 'lO CHANCE UUASAH SELECTION’) 

READ (5,4) nil 

IF(HILNE.XES.ANl).HH.NE.XO) GO 10 74 
IFdiU.EO.XES) lCODE«a 
IE( IIII.E0.X1',S, AND.CC.EO.XFhS) IG0DE*3 
IFdHLEOl.XO.AND.CG.EO.XO) IC0DEa4 * 

IF( ICODE.KO. l.OR. ICODE.EU.O) CO lO 1 
1F( ICOl)E.Eft.4) CO lO 16 
CO IX) 14 

GALL IT.RMI NATION ROUTINE 
CALL FlNllT (0,70) 

S'lOP 

END 


VI* 475© 
VP 476© 
VP 477© 
VP 470© 
VP 479© 
VP 4800 
VP 481© 
VP 482© 
VP 4030 
VP 4040 
VP 4850 
VP 406© 
VP 4H7© 
VP 4800 
VP 409© 
VP 490© 
VP 491© 
VP 492© 
VP 49:il) 
VP 4940 
VP 495© 
VP 49(.0 
VP 4970 
VP 49IM) 
VP 4990 
VP 500© 
VP 5(M« 
VP 502© 
VP 5030 
VP 504© 
VP 505© 
VP 0060 
VP 5070 
VP 300© 
VP 5()9« 
VP 5100 
VP 511© 
VP 5120 
VP 5130 
VP 514© 
VP 5150 
VP 5160 
VP 5170 
VP 5UI0 
VP 5190 
VP 520© 
VP 5210 
VP 522© 
VP 5230 
VP 324© 
VI* 025© 
VP 5260 
VP 527© 
VP 5200 
VP 5290 
VP 5300 
VP 5310 
VP 5320 
VP 0330 
VP 5340 
VP 5350 
VP 536© 
VP 5370 
VP 5380 
VP 539© 
VP 54©0 
VP 841© 
VI* 842© 
VI* 843© 
VP 844© 
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J 

d 






1 


SmiROUTIRB NAPimV (IWTAT» 


HP 

10 

C 


NP 

St 

C 

** 


NT 

80 

c 

** PLOT OlGITiam MAP 


NP 

40 

c 

m* 


NP 

50 

c 

** INPUT t HO TEIWINAL *NO* lU»iPOi«li!E 

** 

NP 

60 

c 

** N8TAT HUraiER OP OTATIONt ON !riy£ 

»« 

in’ 

70 

c 

** AMV TEIWINAL iTATUS AllMY 


NP 

00 

G 

«« 

wik 

NP 

90 

G 

»« OUTPUT i NAP OP UNITED 8TATE8, STATION LOCATIONS. 


NP 

100 

C 

** NENU (SYHIIOL SELECTION) 

«» 

HP 

110 

C 

»» 


NP 

120 

c 

»« calls SUDHOUTINES PiUHE.UNI'ra.RBCT. 


NP 

180 

c 

*« EOUITR.CIRCLE 


NP 

140 

c 

** 


NP 

150 

c 


NP 

160 

c 



M* 

170 

c 

NOTE * NAPDRW MUST BE NDDIPIED TO ACCONODATE OTTIER NAPS - 


NP 

100 

c 

AREAS OP POSSIBLE GRANGE INDICATED IN PHOORAN 


NP 

190 

c 

NAPDRW WILL NOT BE CALLED IP NAP DRAW OPTION NOT SET 


NP 

200 

c 



NP 

210 


DINENSION X( 12) 


NP 

220 


CONWON /DRAW!/ XX( I ) /DRAW2/YY< 1 > /DRAMI/ARAYI 1) 


NP 

230 


DATA ZO/’NOV 


NP 

240 

c 



NP 

250 

c 

DRAV SCREEN NINDOW 


NP 

260 


GALL FIIANE <350.660.270,005) 


Ml* 

270 

c 



NP 

200 

G 

READ DIGITIZED COODINATE^ OP UNITED STATES BORDER PRON UNIT 4 


NI* 

290 

c 

KR I NUNBER OF DIGITIZED NAP RECORDS (GRANGE IP NECESSARY) 


NP 

300 


Kll«77 


NP 

310 


DO 4 J»1.KR 


NP 

320 


READ (4. I) (X(IO,K«l, 12) 


NP 

330 

1 

FORMAT (flX. I2F6.0) 


NP 

340 

G 

PLOT THESE GOORDIN vrVS 


NI‘ 

350 


DO 3 l»l. 11.2 


NP 

360 

G 

SK.TK > TRANSLATION GOWONigNTS (GRANGE IF NECESSARY) 


NP 

370 


SK« 1 1 . 03 


NP 

300 


TK«0.<)2 


Nl* 

350 


X( I)»X( D-8K 


NP 

400 


X( I+I)»X( I+D-TK 


NP 

410 

G 

CONVERT TO VIRTUAL UNiTO FROM CENTIMETERS 


NP 

420 


CALL UNITS (XU) ,X( I+l) .DX.DY) 


NI» 

430 

C 

"PEN UP* FOR FIRST POINT OF FIRST RECORD 


MU* 

440 


IF( I.EO. I.AND.J.EO. 1) GO lO 2 


Ml* 

450 

G 

DIUW TO COORDINATE DX.DY 


Nl> 

460 


CAI.L DRAWA (DK.DV) 


Nl* 

470 


GO TO 3 


NP 

400 

G 

MOVE TO COORD f NATE DX.DY 


NP 

490 

2 

CAIX MOVEA (DK.DV) 


m* 

500 

a 

CONTINUE 


NP 

510 

4 

CONTINUE 


NP 

520 

C 



Ml* 

530 

G 

IWAD DIGITIZED STATION COORDINATES 


Nl* 

540 

G 

DRAW THE STATION NUMUFJIS IN ALPHANUMERIC NODE 


Nl* 

550 


IH) 12 J-l.NSTAT 


MP 

560 


READ (4.5) (X(K) .K«l,2) 


MP 

570 

8 

FORMAT (8X.2F6.3) 


m* 

500 


X( l)»X( D-SK 


NP 

590 


X(2)«X(2)-TK 


MP 

600 


XX( J)«X( 1) 


Nl* 

610 


YV( J)»X(2) 


HP 

620 


CAU, UNITS ( X( 1 ) , X( 2) . DX. DY) 


Ml’ 

630 

G 

DIUW A POINT AT STATION LOCATION 


Nl* 

640 


CALL POINTA (DX.DY) 


MP 

650 

G 

1 

1 

POSITION STA'ITON NUMBERS ON MAP (CHANCE IP NECESSARY) 


Nl* 

660 

IW 
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IF(J.»tt.a> Dy»DY~08. 

MP 

670 


ir< J.Ett.4) DY»DY-15. 

MP 

600 


IFU.EQ.8) 8Y-0Y-20. 

MP 

699 


iru.eo. 1. on, J.EO. 3.011. J. £8.4) D)1>>DX-88. 

NP 

700 


iriJ.Ett.P^ DX«DX4^10. 

MP 

710 


IPU.EO.l^ DX«0X-80. 

Ml* 

720 


IFIJ.E0.4i DX*DX-«-40. 

Ml* 

780 


CALL NOVEA (DX.DY) 

MP 

740 


CALL 8V8TAT ( AIUY) 

MP 

750 


CAU. ANNODE 

Ml* 

760 


IFU.EO. 1> WHITE <6,«) 

Ml* 

770 


IPU.EQ.3) WHITE (6,7) 

HP 

700 


IFU.Ett.O) WHITE 

M* 

790 


IP(J.E0.4) WHITE («.a> 

MP 

000 


iPU.EO.8) Wnim (6,9) 

HP 

010 


lF(J.Ett.6) WHITE (6,11) 

MB’ 

1120 

« 

FORMAT ( • W8*) 

MP 

080 

7 

FORMAT ( • OV ) 

NP 

o-te 

8 

FORMAT < • G8’) 

HP 

use 

9 

FORMAT < • F»’ ) 

MP 

060 

18 

FORMAT ( • CR* ) 

MP 

070 

11 

FORMAT < • HH* ) 

MI* 

000 


CALL RE8TAT (ARAY) 

Ml* 

090 

IS 

CONTINUE 

Ml’ 

900 



MP 

910 


DRAW 8TATI0N SYMBOL NEMU 

MB* 

920 


OPTION TO CHOOSE FROM 8 STATION SYmOLS 

MP 

930 


1. RECTANGLE 

MP 

940 


2. TRIANGLE 

MP 

900 


3. CIRCLE 

MP 

960 


CALL MOVABS <!OCH( 18.8) ,XCN(4.4)> 

Ml* 

970 


CALL SV8TAT (ARAY) 

Ml* 

900 


CALL .\NMODE 

Ml* 
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WRITE (6,13) 

Ml’ 

1000 
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NP 
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CALL HE8TAT (ARAY) 

MP 

1020 


DEFINE NEW WINDOWS FOR SYMBOL DRAWING 

MP 

1030 


CALL SWINIK) (KCM( 12.0) ,KCM(7.») ,KCM(2.8) .RCRIS.O) ) 

Ml* 

1040 


CALL UNITS (7.0,a.0,S,T) 

Ml* 
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CALL VWINDO (0O,8,O.O,T) 

MP 
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CALL FRAFD: (KCM( 12.0) ,)0:M(7.0),KCN(2.a).KCN(2.8)) 
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Ml* 

1000 
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Ml* 
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Ml* 

1100 


CALL EttUrm (3.2, 1.0. 1.0) 

Ml* 
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DRAW CIRCLE 

MP 
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CALL CIRCLE ( 3 . 7, 1 . 0.0. 8) 

Ml* 
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MP 
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** 

SELFATP AND MAP BASELINES 
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US 

c 


INPUT » MAI* 

MAP DRAWING OPTION INDEX 


ns 
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»w 

I CODE 

INITIAL PROGRAM RUN INDEX 

** 

ns 

c 
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NSIAT 

NUMUEll OF STATIONS ON FILE 

** 

BS 
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mm 

BS 
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OUTPUT 1 IN 

NUMBER OF STATIONS CHOSEN 

** 

BS 

c 

♦ ■d 

INO 

NUMBER OF UASELINES CUUSEN 

mm 

IIS 

c 

** 

IS.. IS 

SELECTED BASEL INFil INDEX 

mm 

ns 
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** 

1ST 

SELFX3TED STATION INDEX 

mm 

us 
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*■* 

NEX 

EXl'ElUflENT NUNUEU 

mm 

HS 
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mm 
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** 

OPTIONS 1 STATION 

SYMBOL SELECTION 

mm 

IIS 
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mm 

BS 
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MIK 
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US 
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CALLS SUBROUTINES UN riS.RECT.EUUlTR. CIRCLE 
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INIKGISIl rNU 

niHKNSION IS(I), .IS(I), IST(I) 

COMMON /DRAVI/ XX( l)/DHAVf2/mi)/DRAWa/AIlAYU)/DEXl/IS/DEX2/JS/lMfiX 
la/ IS'1VUEX4/NEX^U«/'TNII, JUST 
Pl»3. 14 

f 1ST AVAII.AIIIJB VI.Bl STATIONS 
CAM. HONK 
CALL AWMOOE 
VHITP, (ft.l) 

CHANGE STATION NAMES IF NECESSARY 
FORMAT ( ’ STATION SELECTION V' 1. WESTFORIIV* 

II.F.YV a. GHF.F.NBAMKV 4. COLOSTONEV 
2 0 . niCUHONDM 


a. ONF.NS VAJ. 
FT. DAVIS V‘ 


STATION AND BASELINE SELECTION 
VIUTE («,2) 

CHANCE NUMBER OF STATIONS AND BASELINES IF NECESSARY 
FORMAT (’ ENTER -^STATIONS, -^BSLNSV STATION 1 2 8 4 8 6V* PUT 0 I 
IF NOT OUSERVINO’.''' MAXIMUM -*'I3SL«S*6’) 

DO a I-I.NSTAT 
ISTf I)«0 

lUCAU 1N,TNU,< 1ST< 1) , 1»1,NSTAT) 

DO f» IJ«1.TNR 
CALL HESTAT (AMY) 

CALL MUVAHS ( I , KCM( 0.0) ) 

CALI, SVSTAT (ARAY) 

CALL ANMORE 
MUTE (6,4) I.l 

FORMAT (' CHOOSE BASELINE #* , lO/* ENTER I J OF OASELIHE*) 

BEAD (8.X0 iS( IJ) , JS( IJ) 

MAP DISPLAY IS SKIPPED IF MAP DRAW OPTION NOT SET 
1F< MAl’.EO.O) CO IX) (> 

CAM, RF.STAT (ARAY) 

CALL MOVABS ( 1 , KCM( 7.2)) 

CAM, SVSTAT (ARAY) 

CAM, ANMODE 

S'lfiHOL SELECTION 
WRITE Mi. 5) IJ 

FORMAT <’ SELECT SYMBOLS '»*,11,* USLINE’/’ PRESS 1- RETURN, MOVE*/’ 


BS 

BS 

BS 

US 

BS 

BS 

BS 

BS 

BS 

BS 

BS 

BS 

BH 

BS 

BS 

BS 

BS 

US 

BS 

BS 

BS 

BS 

US 

IIS 

BS 

US 

BS 

BS 

US 

BS 

BS 

BS 

BS 

US 

BS 

BS 

US 

IIS 

BS 

US 

BS 

BS 

US 

ns 

US 

US 


!• 

20 
ilO 
40 
SO 
()0 
TO 
00 
90 
100 
I 10 
120 
180 
140 
100 
160 
170 

too 

190 

200 

210 

220 

280 

240 

280 

200 

270 

200 

200 

800 

.810 

820 

880 

840 

880 

860 

870 

800 

890 

400 

410 

420 

480 

440 

480 

460 

470 

400 

400 

800 

510 

820 

880 

340 

880 

860 

870 

800 

500 

600 

610 

620 

680 

640 

680 

660 
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ICUnSOR INSIDE SYHhPnESS P*> 

READ (S,«) NON 
CALL RESTAT (ARAY) 

G 

G REDEFINE MENU WINDOW 

CALL SWINDO (KCN( 12.«> .KCM<7.e> .H0M2.8) ,NDM(2.0)> 

CALL UNire <7. ,2. ,8,T) 

CALL VWINDO <«.,S,e.,T) 

C FIND CURSOR LOCATION IN VIRTUAL WINDOW 
CALL VCURSR ( Il.XI.YD 

G SELECT SYimOL ACCORDING TO CURSOR LOCATION 

IF<Xl.Gr.FLOAT<KCH(0.5)> .AND.XI.lt. FLOAT<KCH( 1.0))> I8Y»1 
1F( XI . GT. FLOAi < KCM( 3. 5) ) . AND. XI . LT. FLOAT! KGHC 3. 5) ) ) I8YM>2 
IF( XI . GT. FLOAT! KCH( S . 2) ) . AND. XI . LT. FLOAT! KCH! 6 . 2) ) > ISYN-3 


C REDEFINE MAP WINDOW 

CALL SWINDO ! 350,660, 270,503) 

CALL VWINDO !0. . 1820. ,0. , 1010. > 

C DRAW APPROPRIATE SYMBOL 

C DRAW CHOSEN SYMBOL 

IF! ISYM.Ett. 1) CALL HECT !XX! IS! IJ) > ,YY! IS! IJ> > . 1.0) 

IF! ISYM.E0.2) CALL EaUIlH ( XX! IS! IJ) ) , YY! IS! IJ) > , 1 .0) 
IF! ISYH.Ea.3) CALL CIRCLE ! XX! IS! IJ) ) , YY! IS! IJ) ) , .3) 
IF! ISYN.EO. I) CALL RECT ! XX!.)S! IJ) > , YY! JS! IJ) ) , 1.0) 

IF! ISYM.Ett. 2) CALL EttUITR ! XX! JS! I J) ) , YY! JS! IJ) ) , 1 .0) 
IF! ISYM.Ett. 3) CALL CIRCLE !XX! JS! IJ) ) ,YY! JS! IJ) ) , .0) 

C 

G DRAW BASELINE 

CALL JNITS !XX!1S!IJ)),YY!IS!IJ)),DX,DY) 

CALL MOVEA !DX,DY) 

CALL UNITS !XX<JS! IJ)) ,\Y(JS! IJ)),DX,DY) 

CALL DRAWA !DX,DY) 

1 F! I J . Ett. TNB) CALL MOVABS ! 1 , KCH! 9.0)) 

CALL SVSTAT !ARAY) 

6 CONTINUE 
C 

C DRAW EXPERIMENT NUMBER 

CALL MOVABS ! 1 , KCM! 4 . 4) , DX, DY) 

CALL SVSTAT lAHAY) 

CALI, ANMODE 
WRITE !6,7) 

7 FORMAT !’ ENTER EXPERLTENT #’) 

READ !5,)C) NEX 
IF(M/U*.Ett.O) GO TO 8 

C ALL RESTAT ! ARAY) 

CALL UNITS !9.0, 16.5,DX,DY) 

CALL MOVEA !DX,DY) 

GO TO 9 

8 CALL MOVABS ! 1 , KCM! 2 . 5 > , DX, DY) 

9 CALL SVSTAT ARAY) 

CALI, ANMODE 
WIU'IE !6,10) NEX 

10 FORMAT!’ EXl’EIUMENT . IX, 13) 

C SKIP TO NEW PAGE ON LINE PRirPER IF REPEATING PROGRAM 
IF! ICODE.GT.0) WRITE !7,11) 

11 FORMAT !1H1) 

WRITE (7,12) NEX 

12 FORMAT ! lOX, ’EXPERIMENT , 13) 

G 

CALL MOVABS (KGM(6.;i) ,KCM!3,4) ,DX,DY) 

CALL ANMODE 
WRITE !6,13) 

13 FORMAT !’ PRESS 1 THEN RETURN’) 

READ !5,») MOM 

C 

RETURN 

END 
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BS H40 
BS 1150 
BS 1160 
BS 1170 
BS 1180 
US 1190 
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BS 1240 
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BS 1300 
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READ 

t 

lYEAR. 

. . INITIAL EPOCH OP OBSERVATIONS (1)T> 
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JYEAR. 

.. FINAL EPOCH OP OBSERVATIONS (UT) 
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TO 

too 
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WBITO 
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TOO 

GST AT INlTliU^ EPOCH 
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OUTPUT 
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’TOO 
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180 
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TF 

TOTAL INTERVAL OP OBSERVATIONS 

** 
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140 
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TM 

150 

G 



GALI^ SUBROUTINE CHE8ID.DEGNS 


TO 

'TO 

160 
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170 


C 
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I 

S 
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TN 

PiPLICIT UEAL»0(A“H,L“2i) TN 

rt>NMON /Tim/ INO, IDAY, lYEAR TN 

TN 

HEAD INITIAL EPIKJII IN UNtVEHSAL TINE TN 

WIUTE (6,3) TO 

FOUHAT (• ENTER INITIAL EPOCH IN UNIVERSE > TIME*/’ PORNATi YEAR. NO TO 
INTU.OAY.HOUR.HIN.SEC’/’ PRESS RETURN THEN ENTER FINAL EPOCH IN SIN TO 
2ILAR HANNEll’//’ IF INITIALaFINAL EPOCH IN HIFFEBEHT NONTIIS OR YEAR TN 
»SV EXPRESS FINAL EPOCH IN SANE MONTH OR YEAR AS INITIAL EPOCH’/* TN 
4 E.C. IF INITIAL EPOCH DEC 80, 1979 FINAL EPOCH JAN 1,1900*/’ EN’PE TN 
511 FINAL EPOCH AS DEG 82,1979’//’ ENTER FOR INITIAL EPOCH 0 HOURS U TN 
6T OF INITIAL DAY’) TN 

READ (5,») I YEAR, TOO, IDAY, I HOUR, I MIN, SRC TO 

READ (5,*) J YEAR, .TOO, .IDAY,. IU0UR,JMIN,8ECJ TO 

WRITE (6,8) r 

WHITE (7,3) Th 

FORMAT (//, 15X. ’ INITIAL EPOCH’ ,2X, ’(UT) ’) TO 

WRITE (6,4) TN 

WRITE (7,4) TO 

FORMAT (/lOX, ’ YEAR’ , IX, ’MONTH’ , IX, ’DAY’ , IX, ’HOUR’ ,2X, 'MIN’ ,2X, ’SEC TM 


1 ’ ./) 

WRl'I'E (6,5) lYEAR, IMO, IDAY, IHOUR, IMIN, SEC 
WRITE (7,5) lYEAR, TOO, IDAY, IHOUR, IMIN, SEC 
FORMAT ( lOX, 14 , 4 1 5 , F5 . 1 ) 

WRITE (6,6) 

WRITE (7,6) 

FORMAT (/, 15X, ’FINAL EPOCH’ ,2X, ’ (UT) ’/) 
WHITE (6,5) JYEAH,JM0,JDAY,JII0UR,JM1N,8ECJ 
WRITE (7,5) JYEAll,JMO,JDAY, .IHOUR, JMIN,8ECJ 

IF(JYEAR.NE. IVEAH.OR. JNO.NE. TOO) CO TO 9 

TF I TOTAL INTERVAL OF OBSERVATIONS 


TM 

TO 

TM 

TM 

TM 

TO 

TO 

TM 

TM 

TM 

TO 

TN 

TN 


TF*DFLOAT( JDAY- IDAY) *2*. DO+DFLOAT( J HOUR- IHOUR) +DFLOAT< JNIN- lNlN)/6 TM 


10. DD+( SECJ-SEC) /86O0. DO 

(CALCULATE GIUKENWICH SIDEREAL TIME AT EPOCH TO 
CA1.L CJRESID ( I YEAR. IMO, IDAY, IHOUR, IMIN, SEC, TOO) 

PI2*2.D0»P1 

CHECK FOR NEGATIVE VALUE OF GST 
I F( rilO . CT. “P 12 . AND. TII0 . LT. 0 . DO) TIIO»TO0+P 12 
TII0T»TII0/I5.D0 

CAU, DEGMS ( THOT , P I , I DEG. IMI N , SEC) 

WRITE (6,7) IDEG, IMIN, SEC 

WHITE (7,7) IDEG, IMIN.SEC in 

FORMAT (//lOX, ’GREENWICH SIDEHEAL TIME AT INITIAL EPOCH «’,2I4,2X, TM 
1F6.8) TM 
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420 
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TM 
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READ (S,«) nOH 
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430 TO 11 


TM 
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TM 
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WRITE (6,10) 


TN 

720 

IB 

FORNi\T (//' WRONG INITIAL AND FINAL DATA - RE-ENTER’ ) 


TM 

TN 
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WRITE (6,8) 
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GALL ANMODE 
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** 

»» 
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G 
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** AT INITIAL EPOCH OF OBSERVATIONS 
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** INPUT : lYEAR. .. INITIAL EPOCH OF OBSERVATIONS 

** 

GR 
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** 
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** 
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'** CALLS SUBROirriNE JULIA 
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IMPLICIT REAL!4tO( A-II, L-Z) 


GR 

150 


PIa4.D0*DATAN( l.DO) 


GR 

160 

G 



GR 

17© 


CALL JULIA ( lYEAR. IHO, IDAY, IHOUR, ININ, SEC, NJD) 


GR 
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T« ( MJD-24 15020 . DO) /06525 . DO 


GR 
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CR 

210 
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GR 
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DT»DMOD(BT,24.I)0) 


GR 

230 


CT* ( 0 . 0929DO/&60O . DO) XfT2 


GR 

240 


'niO=A'r+BT+CT 


CR 
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DT=DFL0AT( 1 HOUR) +DFLOAT( IM1N)/6O.D0+SEC/3600.DO 


CR 

260 


rilOs'niO+DT 


GR 

270 


n'(TH0.CE.24. DO) ’ni0»TlI0-24.D0 


GR 
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TU0=rUO* 1,5 . DOK^P 1/ 10© . DO 


CR 
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GR 
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GR 

310 


END 
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INPUT ! lYEAR... INITIAL EPOCH OF OBSERVATIONS 


JL 
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OU'PPUT ! H.I1) JULIAN DATE OF INITIAL EPOCH 

** 

JL 

00 

C 

its* 

** 
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IMPLICIT REAL*0(A-H,L-Z) 
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DIMENSION IMOMIK 12) 
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IF( lYEAR.GT. 1900) ID1S» ID18-1 
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icii«4«( 1YI!;AR/4) 

IF( lYEAR.Ett. iCK.AND. IM.GT.2) ICOD-1 
IF( lYGAH.Ett. 1900) IC01)»0 

HJD«2415e20+< IYEAR-1900>»068.M4-IDIS41NONTH( lN)4I0OD-O.8Dt4lDAY4-H/ 
1 24 . D04^N/ 1 440 . D04^S/06400 . DO 
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SimnOUTINE STATNS (X, Y.Z.XT, YT.ZT.DIflT.PI. 1C0DE,H8TAT> 
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** 
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LONGITUDE OF STATION 

** 

ST 
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II 

HEIGHT OF STATION ABOVE ELLIPSOID 

** 
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** 


IS,JS 

STATION INDICES 

** 

ST 

ST 

100 
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** 
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BASELINE INDEX 
EXPERIMENT NUMBER 
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NEX 

** 
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TND 

lOTAL NUMBER OF BASELINES 

»)|c 

ST 
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NUMBER OF STATIONS ON FILE 
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(CODE 
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ST 
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ST 
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A, FI 

ELLIPSOIDAL PARAMETERS 


170 
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)K* 

ST 

lUO 
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WRITE 

! X,Y,Z 

CARTESIAN STATION COORDINATES (APPROX) 

** 

ST 
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c 



XC,YC,ZG 

COORDINAIE DIFFERENCES 

** 

ST 

200 

c 

!(CD! 


DIST 

BASELINE DISTANCES 

** 

ST 

210 

c 

5(«* 


IDEGI. . . 



ST 

220 
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** 


IDEG2. . . 



ST 

230 
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II 



ST 

240 
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** 


LI,L2,L3 

IS,JS 



ST 

230 
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ST 

ST 

260 

c 

** 




A:): 

270 

c 

** 

OIDTUT 

s XTP.yrp 

BASELINE UNIT VECTOR CONPONENTB 

** 

ST 

280 

c 



Z1T 


«)K 

ST 

290 
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Ht* 


XC,YG,ZC 



ST 

300 

c 



DIST 



ST 

310 
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t* 





ST 

320 

c 





<* 

ST 

330 

c 

** 


CALLS SUBROUTINE RAD 

** 

ST 

340 

c 

** 





ST 

350 

c 


ST 

360 

c 






ST 

370 


IMPLICIT HEAL!JtO(A-II,0-Z) 

INTEUEIl TND 

DIMENSION Xru>. YI’U) , ZT(I), X( I) , YU), ZU), XITU), YTPU), Z 
nv(l), XCU), YOU). ZCU), DISTU), ISU), JSU), ISTU) 

COMMON /cmt/ XTlVCIlDa/YTP/CHDO/ZTP/CIUH/XC/CDDS/YU/CRDt./ZC/'DEXl/ 1 
1S/DEX2/JS/DEX3/ IS'IVDEX4/NEX^DS/TNB, JUST 

DEAD ELLIPSOIDAL PARAMETEIIS PHOH FILE 3 
IFUCODE.NE.O) IlEWIND 3 
DEAD <3,«) A. FI 
WRITE <7,l) 

FORMAT (//, I0X, ’PARAMEIERS OF REFERENCE ELLIPSOID') 

WRITE (7,2) A. FI 

FORMAT </.9X, 'EttUATOIlIAL RADIUS » FI 1 .2,3X, 'METERS' ,/'9X, ' IN’i^ERSE 
1 FLA1TENINC » ' .F7.2/) 

Corn UTE ELLIPSOID FLATrENINC AND ECCENTRICITY 

FLAT* I . DO/F I 

EE»2 . D0»FLAT-FLAT*FLAT 

CALL, EIUSE 
CAI-L HOME 
CALL ANMODE 

READ CHOSEN STATION COORDINATES FROM FILE 3 
WRITE (6,3) 

WRITE (7,3) 

FORMAT (/. 15X. 'AI'PROXIMATE STATION COORDINATES') 

WHITE (6,4) 

WRITE (7,4) 


c 

c 


9 

le 


c 

c 


II 


12 

19 


14 


FOlWAT (/,25X, ’ LATITUDE ' ,4X, ’LOHGITUDE • ,BX, •HEICIIT%/'/,ar 
IX. ’DO* , IX, 'HIK' , IX, ’SEC’ ,6X, ’DEO’ , IX, ’MIX’ , IX, ’SEC’ .9X. 'NE’I1^*> 
INDX>« 

DO 19 |al,NSTAT 

HEAD (9.S) IDUN,L1,L2,L9. IDBOI, Him, SBCI,IDEG2, Hina, SEC2.H 
FORHAT ( I2.9A4. 19, I2.F6.9, 14, 12.F6.9.8X,F19.2> 

IF( I8T( I) .NE. I) GO TO 7 
INDK«INDX4-1 

WHITE (6,6) IDUN.Ll,L2.L8.IDEGI.HINl,8EGI.tDEG2,HIN2,6EGa.H 
WHITE (7.6) IDUH.L1.L2.L9. IDEGl.niNI.SEGl, IDE(X2.HIN2,SEC2,H 
FUHMAT (6X. I2.2X.9A4,4X, 19, IX, 12, lX.F6.a,4X, 19, IX, 12, IX,F6.9,3X,ri 

'convert latitude AJND longitude to RADIANS 
CALL HAD ( IDECI,MIN1,SECI,PUI,PI) 

CALL lUD ( IDEG2,MlNa,SEC2. ALONG, PI) 

COMPUTE RADIUS IN PRIME VERTICAL 
AN«A/1)S0RT( 1 ,D0-EE*D8IN(PHI)*»2) 

COMPUTE CAim:^IAN COORDINATES OF STATIONS 

HA8FJLINE UNIT VECTOR 

Xl'( I ) ■ DC08( PH 1 1 «DCOS< ALONG) 

VT( l)«DG0S(PHI)*D8IN(AL0NG) 

/T( I)»D8IN(PHI) 

X( I)^( AN+H)*X'r( I) 

Y( I ) « ( AN+Il) *YT( I ) 

Z( I ) “ ( AN*( 1 . Ue-EE) +R) *ZT( I ) 

S'lORE CHOSEN BASELINE UNIT VECTORS 
IF( IST( I) .NE. I) GO TO 10 
X1’P( INDX)«X1’( I) 

Y1T( INDX)«YT( I) 

ZTP( INDX)*ZT( I) 

WRITE (6,0) X(1),Y(I),Z(I) 

FORMAT ( I0X, *X « ’ ,F14.9,2X, ’ Y « * ,F14.8,2X, ’Z »’ ,FI4.9,2X) 

WRITE t 7,9) X( 1) ,Y( I) ,Z( 1) 

FORHAT (/,IOX,*X » ’ , F14. »,2X, ’ Y » * ,F14.9,2X, ’Z «’,FI4.8,2X) 
CONTINUE 

COMPUTE COORDINATE DIFFERENCES AND BASELINE LENGTHS 
WRITE (6,11) 

WRITE (7,11) 

FORMAT (/^.6X, BASELINE’ ,6X, ’DELX' , lOX, ’DELY* , lOX, ’DFIJ!:’ ,aX, *DI8TA 
INCE* ,2X, ’iM) •/) 

DO 19 K-l.l'ND 
XC(K)»X( JS(K))“X( I8(K)) 

YC( K) » Y( JvS( K) ) “Y( ISC K) ) 

ZC( K) »Z( JS( K) ) “Z( 1S( K) ) 

D I ST( K) = DSailT( XC( K) **2+ YC( K) *!):2+ZC( K) **2) 

WRITE (6,12) IS(K) ,.IS(K) ,XC(K) ,YC(K) ,ZC(K) ,DIST(K) 

WRITE (7,12) 18(K).JS(K),XC(K),YC(K).ZC(K),DI8T(K) 

FORMAT (7X, I2.2X, I2,4F14.9) 

CONTINUE 

WRI're (6,14) NEX 

FORMAT (///,’ THESE ARE THE STATION COORDINATES FOR EXPERIMENT , 
112/’ PRESS 1 THEN RETURN’) 

IIE,\D (5,*) MOM 

RETURN 

END 


ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

sr 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 

ST 


670 
60S 
690 
700 
710 
720 
790 
740 
750 
760 
770 
700 
790 
BOO 
BIO 
020 
090 
040 
OSO 
060 
070 
000 
090 
900 
910 
920 
990 
940 
950 
960 
970 
900 
990 
1000 
1010 
1020 
1030 
1040 
1050 
1060 
1070 
1000 
109© 
1100 
1110 
1120 
1130 
1140 
1 ISO 
1 160 
1170 
IlUO 
1 190 
1200 
1210 
1220 
1230 
1240 


106 


SUOlWUTiNE 0UA8AR (XD, DECH.ilAH, El, E3. K8. OHG.PI, lOUAi, KEEP. INDEX. 1C Qi! 
tODE, IN. IPFIX.NQUAS.TIW.TF. IN) m 


c 

** 





m 

G 

** 


QUASAR liANDLKR 


QS 

C 

** 




«« 

QS 

G 

** 

INPUT 

1 XTP . VTPj ZTP 
IMO. IDAY. lYEAH 

DAHBLIME UNIT VECTOR 


QS 

C 

** 


DATE 0F START OF OBSERVATIONS 

«N! 

QS 

c 

** 


NQUAS 

IIUNRIUI OF QUASARS ON FILE 

** 

QS 

G 

** 


ICODE 

INITIAL PROGRAM RUN INDEX 


QS 

G 

** 


THO.TF 

GS'IX). INTERVAL OF OBSERVATIONS 


QS 

G 

** 


ill 

TOTAl NUMBER OF STATIONS 

))!» 

OS 

C 

** 


Pl.OMC 

PROGRAM CONSTANTS 


(H<I 

G 

** 





QS 

C 

** 

HEAD 

t IM 

NUMBER OF QUASARS ORSERVED 

** 

(as 

G 

** 


IQUAS 

CHOSEN QUASAR NUNRERS 

** 

PS 

C 

** 


ZNTNAX 

IPFIX 

MAXIMUM ZENITH DISTANCE 

** 

QS 

C 



REFERENCE R.A. QUASAR 

** 

QS 

C 

** 




** 

QS 

G 

in* 

hRITE 

t INDEX 

VlSIlllLITY MATRIX 

** 

(as 

C 

** 


1 

QUASAR KUraiER 

** 

QS 

G 

** 


Ll.L2.La 

QUASAR NAME 

** 

QS 

C 

** 


1 DEC 1 . . . 

QUASAR RIGirr ASCENSION 

** 

QS 

C 

** 


IDEC2. . . 

QUASAR DECLINATION 

** 

QS 

C 

** 


KEEP 

VISIBILITY MA11UX HEADING 

** 

QS 

C 

** 




** 

Q.S 

0 

** 

OUTPUT 

! RA.D 

CHOSEN QUASAR COOIU) INATFJ3 

** 

OS 

c 

** 


E1.E2.E3 

QUASAR UNIT VEC'1X)II COMPONENTS 

t* 

QS 

c 

** 





QS 

c 

** 

OPTIONS 

1 MinUAL VISIBILITY ULTLINER 


(as 

c 

** 





QS 

c 

** 


CALLS SUDROUTINE 

DECMS.IUD 

** 

(as 

c 

** 




■** 

QS 


C 

C 


c 

c 


c 

c 

u 

c 


U8 

IMin.lC IT m:/VL»0( A-H. O-Z) 

INTOIKIl*!! INDEX 

DIMENSION .\TP( 1) . \TT( I) . ZTPU). DECRd). UAUU) . ElU). E2(l), F. 
I'vUl). 11A(1>. D(l). IQUAS(l). KEEini), 5NDEX(6.24) 

COMMON . GHD I / XriV ClUia-' rrP/CRDD/ZTlVT 1 ME/ 1 MO . 1 DAY . I YEAll/MOU 1 /HA/SO 
1U2/D 

CALI. ANMODE 


IF FUOGllAM HE- RUN DO NOT LIST QUASAR FILE 
IF( IC0DE.EU.2.0R. ICODE.F.tt.3) CO TO 0 
WRITE (6. I) 

FORMAT f QUASAR SELECTION') 

WRIIT, ((..2) 

WRI TE (7.2) 

FORJMAT ( nil. 1»X. 'APrUOXIMATE SOURCE COORDINATES'/) 

WHITE (0.3) 

WRITE (7.3) 

FORMAT (20X. 'RICIIT ASCENSION’ .3X. ’ DECLINATION* .//.7X. ,2X, 'NAME' 
1.I3X. 'im ’.'HIN • , 'SEC ,(iX. 'DEC ’.'MIN '.'SEC'/) 

1D:aD QUASAR COORD I NA IT:s 
KM - PACE COUN TER 

KN - NUMHER OF QUASARS DISPLAYED PER PACE OF SCREEN DISPLAY 

KN* 1 

KN**26 

DO 7 1=1. NQUAS 

ID;A1) (3,4) IDIIM.H.L2.L3, IDEGl . MIN 1 , SECl . IDEC2. MIN2.8EC2 
FORMAT ( 12.3A4, 13. I2.F6.0. 14. 12.F6.3) 

IF( I.NE.KN) GO IX) 0 

KM»KM+1 

KW»KN!)!KM 


QK 

QS 

QS 

QS 

QS 

QS 

QS 

QS 

QH 

QS 

QS 

QS 

QS 

QS 

QS 

QS 

QS 

QS 

QS 

QS 

QS 

QS 

QS 

QS 

QS 

QS 

QS 

QS 

QS 

QS 

QS 


107 


le 
26 
36 
46 
56 
66 
76 
00 
66 
166 
1 16 
126 
136 
140 
156 
166 
176 
106 
166 
266 
216 
226 
236 
246 
250 
260 
276 
200 
266 
306 

:no 

326 

336 

340 

356 

3«6 

37 ' 

: .» 
366 
466 
416 
426 
436 
446 
4.16 

466 

476 

406 

460 

.166 

SI6 

526 

.136 

146 

116 

.166 

176 

106 

166 

606 

616 

620 

(.36 

646 

656 

666 


ftft fs ana 



WHITE (6,42) 

OS 

670 


RFAD (3,w) MON 

OS 

600 


CALL EHA8E 

US 

690 


CALI, IIOPIK 

OS 

700 


CALI. ANMODE 


710 

8 

WlUIE (6.6) 1,L1,L2,LU, li)i!»;i,niMl,BECl,IDECai,HlN2,BEC2 

US 

720 


WHITE (7,6) 1. 1.1. 1.2. I„T. IDF.ni, MINI. BF.OI, IDEfi2.Nllf2,NF.G2 

OS 

730 

6 

FOnHAT (6X. I2.2X,3A4,4X, 10, IX, 12. 1X,F6.0,4X. 10, IX. 12, lX,F6.8.aX,Fl 

US 

T40 


10.2) 

OS 

730 


CONVERT HA AND DFX! TO RADIAIfH 

Util 

760 


CALL RAD (IDEGl.NINl.BEGl.lUUUD.l'l) 

US 

770 


HAIK |)«UAIl( l)*10.DO 

OPI 

700 


CAM, RAD ( IDF.C2.NIM2.SBG2,DECR( 1) .PI) 

OS 

790 


CD«DCOS(DECiU D) 

US 

000 


WU<I)SIN(UKCR( D) 

OS 

010 


CA«DC08(RAR( D) 

OS 

020 


fc!A«USlN(HAR( D) 

US 

030 


CONPUI'E CONPOHENTR OF OIIASAH UlflT VECTOR 

OS 

040 


El( I)«CDiPCA 

US 

050 


E2( l)*CU*SA 

US 

nr>o 


E!l( 1)»SD 

OS 

070 

7 

CONTINUE 

US 

uuo 



OS 

090 


WRITE (6,42) 

OS 

900 


READ (3,!K) non 

US 

910 

B 

CALL ElUHE 

OS 

920 


CALL HOME 

OS 

90© 


CALL ANMODE 

US 

940 



OS 

900 


DO 9 J» 1 . 12 

US 

960 


UA( J)=O.DO 

US 

970 


IX .DxO.DO 

OS 

900 

9 

CONTINUE 

US 

990 



OS 

1000 

10 

WRITE (6.11) 

08 

1010 


CHANCE FORMAT IF MAXinUN "'i 12) OF ttUABARB INCREAHED 

US 

1020 

11 

FORMAT (//',’ 1M)W J1ANV OUASARS DO VOW WISH TO OBSERVE* ’ NAXtHUN 

OS 

1030 


INUMRER IS 12*) 

08 

1040 


READ (3,!*) IM 

OS 

1030 


CHANCE IF MAXIinm #( 12) OF ailASARS INCREASED 

OS 

1060 


IFUM.GT. 12) GO TO ID 

US 

1070 



OS 

tOBO 


OPTION TO CHOOSE OUASAHS WITHOUT VISIBILITY 0L«1LIHER 

OS 

1090 


IF IFY»0 SKIP VISIDlLm OUl'LlNER 

f(is 

1100 


IFY=(* 

OS 

1 1 10 


WRITE (6,12) 

OS 

1120 

12 

FORMAT (5K, 'DO YOU WISH lO CHOOSE UUASARS DEFORM VISIBILITY OUaiN 

OS 

1130 


IF.U ?' ) 

OS 

1140 


READ (5,13) GO 

US 

1150 

10 

FORMAT (A4) 

US 

1 160 


IFICC.EO.KO) CO TO IB 

US 

1 170 

14 

WRITE (6,15) 

US 

UUO 

13 

FORflAT (/5X, ’ENTER CUOSEN ttUASAIlS V3X, ’E.G. 1 3 3 7’) 

OS 

1 190 


READ (5.*) ( IOUAS(K),K<*l, IN) 

US 

1200 


DO 16 Ksl.lN 

US 

1210 


CHANCE IF NUMDEH OF OUASAHS (47) IN FILE CHANCED 

OS 

1220 


IF( lOUAS(K) .GT.47) GO TO 14 

OS 

1230 

16 

CONTINUE 

OS 

1240 



OS 

1230 


INDX : CUOSEN QUASARS INDEX 

US 

1260 


INDX’O 

OS 

1270 


IFY« 1 

OS 

12B0 


KL«0 

US 

1290 

17 

KL»KL+I 

OS 

1300 


INl)aIOUAS(KU 

OS 

1310 


IF(KL.LE. IH) GO TO 01 

OS 

1320 


CO TO 33 

OS 

1330 



OS 

1340 


108 


18 

CALL F41A8E 

OS 

I3B0 


CALL HOUR 

yg 

1369 


CALL ANItMUR 

m 

1070 

C 


os 

I3B0 

C 

MUTUAL VISIBILITY OUTLIIfEB 

QS 

1390 

C 

UUAV TTTLL KOU VISIUILITV UUTLINEII 
WIUTR (6,19) 

os 

1400 


QS 

1410 

19 

FORMAT ( ' MUTUAL VISIBILITY OUTLINEHV/, • IHPUT NAXIHUn ZEHITH BIS 

US 

1420 


ITANCE’ ) 

QS 

1400 


WITR (7.20> 

QS 

1440 

20 

FORMAT < lili. IOX. 'MUTUAL VISIBILITY OUTLINEH'/) 

QS 

1400 


WIUTK <7.2B> IIM). IDAV, IVBAR 

QS 

1460 

C 


QS 

147© 

C 

INPUT flAKIMUH OBSERVABLE ^ENITU BISTANCE 

QS 

l4tM) 


HEAR («.») /NTMAX 

QS 

1490 


WRITE (7,21) ZNTMAX 

QS 

1000 

21 

FORMAT (/>',5X, 'MAXIMUM ZEN ITU BISTANCE »',F6.2/) 

QS 

1010 

fi 


QS 

1520 

C 

COMPUT'E VISIBILITY MATRIX OUASAIl BY QUASAR 

US 

1030 


1)0 22 IJ’^1,24 

QS 

1540 

22 

KRRP( U)»I.I“I 

QS 

1550 


TO“O.BO 

US 

1060 


INDN) 

QS 

1570 


INBX«0 

QS 

15U0 


1)0 04 1«1.NQUAS 

QK 

1590 


INI)« INBHI 

QS 

1600 


DO 24 K> 1 , 24 

QS 

1610 


T^TO+iK-l) 

QS 

1620 

C 

TRANSFORM OUASAR imiT VECTOR TO EARTH FIXED SYSTEM AT EPOCH T 

QS 

1630 


ANGLE® OMG>l* ( T-TO ) +TUO 

US 

1640 


CA=UCOS( ANCl.E) 

OS 

1600 


SA«nSIN( ANGLE) 

QS 

1660 


QI®CA«<EI( l)+SAit«E2( I) 

QS 

1670 


02=“SAtEl( l)+CA*F,2< 1) 

QS 

I6R0 


O0*E0( I) 

QS 

1690 

C 

CALCULATE ZENITH BISTANCE 

QS 

1700 


1)0 20 ,I» 1 , IN 

QS 

1710 


ARGU*XTP( J) *Q1+YTP( J) *Q2+Z1T( J)*QO 

QS 

1720 


1 F( AUCU . LT’. - 1 . JIO) AUGUs - 1 . 1)0 

QS 

1730 


T F ( AHCU . GT . 1 . DO ) ARGU* 1 . DO 

QS 

1740 


ZNT*l)ARCOS(ARGU) 

QS 

1700 


ZNT= l)ARS( ’/NT) * IMO . I)0/P 1 

QS 

1760 

C 

IF ZNT.GT. ZNTMAX CAUSE TWO STARS TO BE PLACED 
IN APPROPRIATE LOCATION OF VISIBILITY MATRIX 
1 F( ZNT. <; r. ZN TMAX) ’/NT» 1 000 . 1)0 

QS 

1770 

C 

QS 

1700 


QS 

1790 

C 

FILL VISIBlLI'rV MATRIX 

QS 

lUOO 


IN0EX( J , K) “'/NT 

QS 

litIO 

22 

CONTINUE 

QS 

1820 

24 

CONTINUE 

QS 

KK'O 

€ 

DRAW TTILF, FOR VlSiBlU’I'V MATRIX 

QS 

1840 


wn 1 I F. ( (. , 25 ) I MO , I DAY. I YEAR 

QS 

1850 

2S 

FORMAT {^^/IX. 12,’/ ',10,'/", 15, IX, 'ZEN mi BISTANGES («* DENOTES NO 

QS 

UI60 


INVISIBII.ITV) *) 

QS 

1870 


10UTE ( 6 . 26) ( KEEP( ID , IL« 1 , 24) 

QS 

1800 


VRITE (7,20) (KEEP( ID , IL«1,24) 

QS 

1890 

26 

l ORMAT (//.2X, 'OS ST* ,24( IX. 12) ) 

QS 

1900 


DO 20 IK“I,IN 

QS 

1910 


WIUTE (6,27) IND, 1K,< 1NDEX< IK.KI) ,KI»I,24) 

QS 

1920 


WRITE (7.27) INI). IK.( INDEX( IK,KI).KI»I,24) 

QS 

1930 

27 

FORMAT ( IX, 10, IX, 12, 24( IX. 12)) 

QS 

1940 

2U 

CONTINUE 

QS 

1950 


WRITE (7,29) 

QS 

1960 

29 

FORMAT (/) 

QS 

1970 

C 

CHOOSE QUASAR IF APPROPRIATE 

QS 

1980 


wnrrF. (6,oo) 

QS 

1990 

ao 

lORMAT (/’ DO YOU WISH IV OBSERVE THIS QUASAR’/' ANSWER YES OR HO’ 

QS 

2000 


I/) 

QS 

2010 


HEAD (S. 13) GC 

US 

2020 


109 


ftO 


C 

31 


32 


33 


34 

C 

39 


36 


37 

3U 

r. 

39 

43 


41 

42 
C 


(UJASAnS OUT OF M2,' QUASARS') 


iriOG.EQ.KO) CO TO 33 
AUl) OllOSliN UUASAll 
IN3K-IN3X+I 
RA( INDX)»RA1U INI)) 

D( INIWxOI^OiU INI)) 
lOSJASt INDN)«IND 
IF( IFY.EQ. I) CO TO 17 
WIHTK <6,32) INDX, IN 
FOIINAT ( ' YOU HAVE CHOSEN’ . 12. ’ 

VIUTE (6,42) 

KKAI) (5,*) MOM 

IF( INDK.EQ. IM) GO TO 39 

CALL EtlASE 

CAM. HOME 

CALL ANMODE 

CONTINUE 

CALL ERASE 
CALL HOME 
CALL ANMODE 

I, I ST CHOSEN OUASAHS 
WRITE (6,36) 

WHITE (7,36) _ 

FORMAT < 1 HI. //,20X, 'THESE ARE THE SOURCES CHOSEN' /.29X. ’RIGHT ASCE 
INSION' ,3X, 'DECLINATION', //,20X,’^ ’,4X, ’UR ','MIN 
2’, 'MIN ’,'SE(:'/) 

DO 30 l» ( , INDX 
RQ>^RA( D/19.D0 

CALL DECrW ( IMI, I* I , I DEC I , Ml N I , SEC I ) 

CALL DF.CMS ( D( I > . 1> I , IDEG2 . MIN2 , SEC2) 

WRITE (6,37) I, IOUAS( I) , 1UEC1.N1N1,SECI, IUEG2.NIN2,SEC2 
WHITE (7,37) I, 1«UAS( I), IDECI ,MINI .SEC1 , IDEC2, MIN2.8FX2 
FORflAT USX, 12,'. ’ , I2,3X,2( 13, IX, 13, 1X.F7.3.4X) ) 

CONTINUE 


WRITE (6,40) 

FOllflAT </, I9X, 'ENTER THE »U*:FRIU''.NCE SOURCE’/, I9X, 'USE SEQimNTIAL 
IMRER (LEFT- MOST COLUm ABOVE)') 

READ (9,'t‘) ITFIX 

CHANCE IF MAXIMUM .•'I 12) OF OUASARS INCREASED 
IF( IPFIX.CT. 13) GO TO 39 
WRITE (7,41) ll'FlX 

FORMAT (/2(X, 'S0IIIW:E' , 14.' IS THE REPERli.NCE SOURCE') 

FORTLAT (//• PRESS 1 - THEN RETURN') 

RETURN 

END 


08 

(Its 

OS 

OS 

OS 

OS 

OS 

OS 

Oii< 

US 

OS 

OS 

OH 

OS 

OS 

OS 

OS 

OS 

OS 

OS 

OH 

OS 

OS 

OS 

OS 

OS 

OS 

OS 

OH 

OS 

OS 

OS 

OS 

OS 

OK 


OH 
NU OS 
OS 
OH 
OS 


Of« 

OS 

OS 

OH 

OS 

OS 


as30 

2040 

2090 

2060 

2070 

2080 

2090 

2100 

2110 

2I'20 

2130 

2140 

2190 

2160 

2170 

21(10 

2190 

2200 

2210 

2220 

2230 

2240 

2250 

2260 

2270 

22U0 

2290 

2300 

2310 

2320 

2330 

2340 

2350 

2360 

2370 

2380 

2390 

2400 

2410 

2420 

2430 

2440 

2490 

2460 

2470 

24(10 

2490 

2900 



1X0 


i 


« 

G 

C 

G 

G 

G 

C 

C 

G 

G 

G 

G 

G 

G 

G 

G 

C 

C 

C 

G 

C 

G 

G 

C 

G 

G 

C 

G 

C 

G 

G 

C 

G 

G 


f 


I 

I 

t 

1 


C 

I 


2 


3 

4 
C 

n 

6 

G 


simRoirriifE sinult (X£nN.xiLE,XEs,eiiAi»,pnx.rar,P3iz,oi».Tr,im.iriij; 

1, IGOWE.GOWV) 


K«*|i 




** 

*4 


SGimDULE SIMULATOR 

** 






4>iH 

INPUT t 

OMO 

EARTH ROI'ATION lUTE 

** 

!):*( 


EIUD 

MEAN EARTH RADlin? 

Wik 



THO,TF 

GS10. INTERVAL OF OBSERVATIONS 




IGODE 

INITIAL PROGRAM HUN INDEX 

** 

** 


TNB 

TOI'AI. NUMBER OF BASELINES 

** 

!(flK 


ISTT 

NUMBER OF EARHl OHIENTATION STEPS 

** 

Di* 


STSZ 

END OF STEPS KFOGH 

** 



XG.VC.ZC 

RASEUNE COORDINATE DIFFERENCES 

** 

** 


RA.D 

OUASAR GOORUINATl!» 


If:* 


GONV 

GONVFJIT UNIVERSAL ’YO SIDEREAL TIME 


** 




** 

If* 

READ s 

PMX,PMY,PMZ 

Al'PUOK. VALUtaS EAll’lU ORIENTATION 

** 

** 


DT 

TIME INTERVAL BETWEEN OBSERVATIONS 


*«< 


IFILE 

OBSERVATION STORAGE FILE NUMBERS 

** 

** 




** 

** 

WRITE ! 

PMX, PMY.PMZ 


** 





** 

** 

OUTPUT 1 

SIMULATED OBSERVATIONS ON FILES 10-10 

** 

** 


L 

BASELINE NUMBER 


** 


K 

OUASAB NUrWER 

** 

** 


(HOUR, IMIN 

KPOGH OF OBSERVATION 

** 

** 


DS.FRNG 

OBSERVED DELAY Q DELAY RATE 

m 

*!K 




** 

** 

OPTIONS s 

SIMULATE DELAY RATE OBSERVABLE)^ 

Mttff 

** 


READ OBSERVATION SCHEDULE FROM FILE 9 

** 

** 


1U:AD OBSERVA'ITON SGIIEDULIS FIM)N TEHMINAl, 

** 

** 


INPUT TIMES OP OBSERVATION ( DT NOT CONSTANT) 

** 





** 




ira*LIO£t KE/VLH<0<A-U,0-E) 

iNTKGKit Trm „ 

niHRfWION IHI). lUCI), XC(1). YCU). ZC( 1) , PMX(l), pmf< I) , 
1, IFILEU), S’i'SZU) 

GOMflON /V.im/ XG/Gnnn/VG/GRnf»/ZG/ROtJ I -^RA/SOUa/n/BS/TWB , JUST 


PHZ( 1) 


WUITE (0,1) 

rOHNAT (/,• READ tN APPROX VALUER FOR EARTH ORIENTATION '/^ THERE 
1 SHOULD BE a>l« 'STEPS) OF PAMHETERSV FORMAT i TWO POLAR MOTION 
2G0M1‘0N)«;H’1-S 1R MEIERSV EAJl'lH RO'I'A’IION IN KEGONDR OF TIME’) 

READ ( 5 . *) ( PMX( .1) . J» I . ISTT) , ( PMV( J) . J» I , ISTT) . ( PMZ( J) , J» I . I8TT) 
imi TE < 7 , 2) ( PfIX( J ) , J« 1 , ISIT) 

FORMAT (//9X, 'APPROXIMATE VALUES FOR EARTH ORIENTATION' //20X, ’STEP 
1 1 ' . 3 X, ' STF.P2 ’ . BX, • STEPS ' , BX. ' STEP4 ’ / 1 IX, ' PMX* , IX, ’ ( METERS) * , 3X, 4( F 
20.U.2X)) 

Wll I I’E (7,0) ( PNV( .1 ) , ,l» 1 , I sri’) 

FORMAT < /- 1 1 X, ' PMY' , 1 X, ' < mri ERS ) ' , OX, 4( FO . 3 , 2X) ) 

WRITE (7,4) (PMZ(J) ,J»i, ISri) 

FORMAT ( / 1 1 X. ' PM7, ' , 1 X. ' ( SEGONDS) ' . «X, 4( F8 . 8 . 2X) ) 


WRITE (6,0) 

FORMAT (//•' ENTER TIMf- INTERVAL BETWEEN OBSERVATIONS (IN MINUTES)' 
I/' INPUT 0 IF TIME INTLRVALS NOT REGULAR’) 

READ (0,*) DT 
I)T*DT/60.D0 
WRITE (6,6) 

FORMAT (//• PRESS I - THEN RETURN') 

READ (B,*) MOM 

IFdGODE.EO.O) GO TO 7 
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SH 

SM 
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SM 

SH 


It 

at 

at 

40 

(ft 

60 

70 

Rt 

00 

100 

iio 

120 

too 

140 

too 

Hit 
170 
I IK) 
190 
200 
‘210 
220 
2U0 
240 
250 
260 
270 
200 
290 
300 
3(0 
320 
330 
340 
330 
360 
370 
3H0 
390 
400 
410 


420 

430 

440 

400 

460 

470 

4)10 

490 

000 

510 

320 
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540 
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070 

000 

090 

600 


610 

620 

630 

640 
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ntn ftfton ocin a naa nnaa 


C 


7 


G 

C 

U 

9 

10 



la 


10 

14 


in 


10 

17 


10 

c 


lURVIHD IIIPOT riLM 
MW1W0 9 
CAI4, KIIAHK 
CALL ROME 


ir RE-RDIfNIItO PnOORAN 


imiTP <0,R) 

CHANGE FORMAT IF NAXIHON H6) OF DASBLINEi 1NGIIBA8ED 
FUlOIAT (/' CUOOHE ONE FILE PEil OA0BL1NEV AVAILABLE PILE NOHMIllfl 
It 10- in RTART WITH 10 ') 

V11ITK (0, 10) 

FORMAT </• ENTER OUTPUT FII.PJI') 

RF.AI) (5,*) ( IFILE(J) .J*t.TNB) 

DO 11 J^liHIO 

CIIANI.T, IF MAX i NUN -flO) OF RARF.I.INPil INCIWABED 
IF( IFILE(J) .LT. 10. OR. IFILECJl.GT. 15) CO TO 9 
CONTINUE 


IF( ICODE.EC.O) 
DO la l>l,TNU 
NUM'iriLEt I) 
REWIND HUM 
CONTINUE 


CO TO 10 


OPTION TO ENTER ODBERVATIONB AT TERMINAL OR FROM FILE 

1. I FLAG- I t FROM TF.BHINAL 

2. IFLAC*2 t FROM FILE 9 
WRITE (0, 14) 

FORMAT (/• IF VOU WISH TO ENTER DATA AT TERMINAL. INPUT TERMV IF 
ITOU RAVE STORED OBSERVATION SCIIEDULE DATA OH FILE. INPUT FILE’) 
READ CJ,17> CG 

IFIGG.NE.XERW. AND.CO.NE.XILE) GO TO 10 
I n GU . Ett . KICRN) I FLAG* 1 
IFICG.EU.XILE) IFlJlOa 

OPTION t IPASSoO SIMULATE OUSERVATlONi^ EVERY DT MINUTES 

1 PASS" I SIMULATE ODSERVATIONS UNEVEN INTERVALS 

IPA88«0 
WRITE (0,10) 

FORMAT ( • ARK 0RBERVATI0N8 AT UNEVEN INTERVALS OF TIME* ) 

READ (5.17) GG 

LAST OUb’ERVATION SHOULD UE GREATER TRAN TF IF I PASS* I 
IF(GG.EQ.XF.S) 1PASS*I 

OPTION t IFRNOO - SIMULATE DELAYS ONI.V 

IFRNG-1 - SIMULATE DELAY RATES 100 

UTING«0 
FHNfi-O.DO 
WRITE (0, 10) 

FORMAT (/' DO YOU WlSU lO SIMULATE DELAY 1UTE8’) 

READ (n. 17) F 
FORMAT (A4) 

U’U' FO.XKS) 1FRNC»1 

r t TIME COUNTER 

IJ ! H'l’KP COUNI’KII 

1. I BASELINE COUNTER 

T*O.DO 

L)»0 

L«0 

CALL ERASE 
CAI.I, IlONF. 

CALL ANMODE 

UPDATE EARTH ORIENTATION STEP 
IJ* iJ+1 

( IIANCE UN ns OF FJUITH ORIENTATION VAT.UF4I FOR THE I.I TH STEP 
PMI«PMX( M)/ERAD 
Pfl2»PMY( IJ)/ERAD 


•H 

SM 

SM 

SH 

SM 

SN 

SM 

SM 

SM 

SM 

SM 

SM 

SM 

SN 

SM 

HN 

SN 

SM 

SM 

SN 

SM 

SH 

SM 

SN 

SH 

SH 

SN 

SH 

SM 

SM 

SN 

SM 

SM 


070 

000 

090 

700 

710 

720 

700 

740 

750 

700 

770 

700 

790 

000 

010 

020 

000 

040 

ono 

000 

070 

000 

090 

900 

910 

920 

900 

940 

950 

900 

970 

900 

990 


SH 1000 
SM 1010 
SH 1020 
SN 1050 
SN 1040 
SN lono 
SM 1000 
SM 1070 
SM 1000 
SN 1090 
SN 1100 

SM me 

SM 1120 
SN 1130 
SM 1140 
SN lino 
SN 1100 
SM 1170 

SN lino 

SM 1190 
SN 1200 
SM 1310 
SM 1220 
SN 1250 
SN 1240 
SM 1250 
SN 1200 
SM 1270 
SM 1200 
SN 1290 
SM 1500 
SN 1510 
SM 1520 
SM 1550 
SM 1540 
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m>PlfZ( IJ)«15.D0/'30«SMS.M> VN 

IJ) gH 

ir(T.CJT.TP) CO TO ST «M 

irUl'M^.EU* (.AKO. IJ.CT. 1) 00 TO SI ON 

INCIlKAtiK HAOKLINK COUNTER NH 

I.-I.+ I SN 

im..CT.TNU} CO TO 0« SH 

fP(i..cT. I) CO TO oa SH 

ircT.EQ.O.DO.On. IFLAG.E0.3) 00 TO 33 ON 

CALL KIlAOli; SH 

CALI. HOMF. 8H 

CALL AHMODE ON 

r*rf 

vniTP. pnr.vioi»i orneuvation on ocreen sn 

WUITE C6,2i> K, lUOUn.NfN ON 

KOliNAT (’ TIIK PIIKVIOUH ORORHVATION WAN (V OUANAR >M8,3X. ’ HOOR HH 

I «• . I8.2X, • MIN . IS//) ON 

ON 

|P( IPAOO.KO. I) CO TO a* 8N 

IP< IFLAC.Ea.2) GO TO 29 ON 

ON 

r.NTP.R OUARAR NUMBER SN 

WHITE (6,24) OM 

l OlUUT (OX. •CHOOSE NEXT OBSERVATION* /5X, ’ INPUT #<UJASAR* > SM 

READ (O.NO K SN 

CHANCE IF NAXINUN ^( 12) OF UUASAIIO IMCREAOEII OM 

lIMK.CT. 12) CO TO 20 SN 

CO TO 01 SN 

UEAB C),*) K. SN 

CO TO 01 SN 

OH 

1F( 1FLAC.E0.2) CO TO 20 HN 

ENTER (lOASAR AND EPOCH OP OBSERVATION SN 

WRITE (6,20) OM 

FOllNAT (9K, ’CHOOSE OUSEUVATK'N SCHEDULE VOX, * INPUT #OUASAR ROHR N SN 
IINinT.*) SN 

lU-IAD (5,!I<) K, IUOUH.MIN ON 

CHANCE IF NAXINUN ■»( 12) OF OUASABS INCREASED SN 

IF(K.CT. 12) CO TO 2T SN 

CO TO 00 SN 

IF(T.CE.TF) CO TO 07 SN 

READ (9.*) K. 1H0UR,MIN SN 

T« DFLOAT( 1 HOUR) + Mi N/60 . DO SN 

IF(T.CT.ST) CO TO IB 8M 

SN 

CALCULA'IE I'RICONONETHIC HF^NDERS OF OBSERVATIONS SN 

cn«ncos{ D(K) ) sn 

SD=USIN(U(K)) SN 

V I = ONC*'l’+ rilO+PNO*COHV SN 

W-YI-RAIK) SN 

CKI'»UCOS( V2) SN 

SK1'»!)SIN( V2) SN 

SN 

DELAY OUSEUViUlLE SINULATION SN 

DR« - XCn .) *( CD*CKP+PMI #SD) +YC( L) *( Cr^S?KP+PM2K!SD) -ZC< L) *< SD-PHDKCIMtG SN 
IKI*-Pm*i-0D«t8KP) ON 

IF( IFRNC.NE. I) CO 'TO 00 SN 

SN 

DELAY RATE OUSERVADLE SINULATJON SN 

FRN(; IS DERIVATIVE OF DELAY W.H.T. TINE ( HETERS/HOim) SN 

FRNC»ONC*CD*i:( XC( L) *SKP+YC( L) n<CKP"ZC( L) *< PH1>«SK1*-PN2D*CKP) ) SN 

SN 

.«'mnE SIMULATED OBSERVATIONS ON UNITS 10-14 SN 

IF( 1PASS.F.Q. i) CO TO 34 SN 

IHOUR* IDINT( T) SN 

TMIN»T-illOUn SN 

TN 1 N» TM I N t.60 . DO+0 , 0000 IDO SN 

NlN'i IDINT( TMIH) SN 
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ir<Nin.R£.6e) go to 94 

em 2000 


P1IW«MIW“60 

RN 2040 


iHoim«inouR+t 

8H 20S0 

94 

9im«IFILEiL> 

ON 2060 


(HUM. 90) L,K, lUOUn.mN.DG.riUfC 

RN 2070 

99 

FOtUI4T (4IS.3F20.I9) 

SN 2000 


r.0 TO 20 

9N 2000 

96 

l,«0 

RN 2100 


IF( IPASS.EO. 1) GO TO 20 

8N 2110 


T»'|+1»T 

RN 2120 

C 

CIIEGK FOR STEP UPRATF. 

RN 2100 


IF(T.GE.ST) GO TO 10 

OH 2140 


GO TO 10 

RN 2ino 

C 


SH 2160 

97 

RETURM 

RN 2170 


EHU 

RN 2100 
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C 

V 
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G 
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(;UU1U>UTIKE FLAGlj (JK.) 

itHKilo 

l£Xi*EIUNENT FLAG UANDLliK 
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2 


G 

G 
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4 


G 

G 


C 

G 


5 

6 
7 
» 
9 


G 


H:iK 

** 

>N< 

*>K 

it!* 

iKiX 

iK* 


HEAD s FLAG1« It DELAY Uii ONLY OUSEIIVAULE 

FI.ACf«2tOnSEIWAI»LP.S ARE DEI.AV AND DP.T.AY RATE 
FLAGl-3t DELAY RATE IS ONLY OBSERVABLE 

FI.AGa* » I MWl.TI-RASEI.JNE GONF IGURATION 

FLAG2« 2 t ESTIMATE KBl GONPOHENT OF POLAR NOTION 

FLAGa»*Jt ESTIMATE E'PA GOfU'ONBNT OF POLAR MOTION 

FLAGO« It COVARIANCE ANALYSIS ONLY 
FLAGtl»aiGOMl*LETE LEAST SOUARES ESTIMATION 

FLAG4* It ESTIMATE ALL i UUVMETERS 
FLAG4* 2 1 DELETE GI.OGK PARAMETERS 


** 

** 

UiNi 

!«* 

** 

IKIti 

*i|i 

** 


le 


WRITE t FLAG »1ESSAGES 
OUTPUT t JK - ONE BASELINE CASE INDEK 

I NTEGER FLAGl , FLAG2 , FLAGJI , FLAC4 
COMMON /FLG/ FLAGl ,FLAG2,FLAG0,FLAG4 
imiTF, ((.,2) 

FORMAT (5X. 'CHOOSE EXI'EIUMENT FLAGS’/’ FLAG1«! 1 i DELAY IS ONLY 01H$ 
lERVADLE’/’ Fl,AG1=2 t DEI, AY RATE OBSERVABLE INCLUDED’/’ FLAGI*3 t D 
2ELAY RATE IS ONLY OBSERVABLE’//’ FLAG2» 1 t MULT I- BASELINE EXPERIME 
UNT’/’ FLAG2=2 ONLY ETA GOm’ONENT OF POLAR MOTION’/’ PLAC2*3 t ON 
4LY KSI GOMPONEiIT OF POI.AR MOTION’//’ FLAGS* I : COVARIANCE ANALYSIS 
N ONLY’/’ FLAG3=2 t COMPLETE LEAST SOUAIUSS SOLUTION’//’ FLAG4* 1 I A 
6LL I’AIIAMETEBS’/’ FLAG4=2 t NO CLOCK PARAMETERS’//’ INPUT FLAGl 
7 . FLAG2 , FI.AG3 . FLAG4 ’ ) 

lUOAJ) PROGRAM FIvAGS 

READ («.*) FLAG1,FLAG2,FLAG3,FLAG4 

1F(FL;-G1.GT.3.0ILFLAG2.GT.3.01LFLAG3.GT.2.0R.FLAG4.GT.2) GO TO 1 
WHITE (7,3) 

FORMAT (//, 12X, 'PROGRAM FLAGS’ ) 

WUl'l'E (7,4) FLAGl , FLAG2 , FLA(JU , FLAG4 

FORMAT (/,9X. ’FLAG! » ’ , I2/9X, ’ FLAG2 » ’ , I2/9X, ’ FLAG3 » ’,I2/9X, 'P 
1LAG4 * ’ , 12//) 

IF ONE BASELINE OMIT ORE POi 1R MOTION PARAMETER! JK* I ) 

JK*0 

1 F( FLAG2 . EO. 2 . OR. FLAG2 . Eft. 3 ) .) K= 1 

Wlin'E FLAG Mi:SSAGES 
IF(KLAGI.EO. 1) WRITE (7,5) 

IF(FLAGl.Ea.2) WRITE (7,0) 

1F(FLAG1.E(A.3) WHITE (7,7) 

IF(FLAG3.EQ. I) WRITE (7.8) 

IF(FLAG3.Ea.2) WRITE (7,9) 

FORMAT (9X, ’ANALYSIS INCLUDES ONLY THE TIME DELAY OBSERVABLE’) 
FORMAT (9X. ’ANAI.YSIS INCLUDES TIME DELAY8T1ME DELAY RATE’) 

FORMAT (9X, ’ANALYSIS INCUDES ONLY DELAY HATE OUSEllVAHLE ’ ) 

FORMAT (/,9X, ’COVARIANCE ANALYSIS ONLY’) 

FORMAT (/.9X, 'COMPLETE LEAST SOUARES SOLUTION’) 

WHITE <0,10) 

FORMAT (//’ PIUISS 1 - THEN RETURN’) 

READ (0,«> NON 

RF/niRN 

END 
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190 
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000 
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040 

050 

000 

070 

680 
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OUSERVATION WEIGHTING HANDLER 


wr 

** 





WT 


INPUT 

> 

TNB TOTAL NUMBER OF BASELINES 


WT 

HiiH 



FI FLAG FOR OBSERVABLES 

** 

wr 


READ 

I 

SIGI,SIG2 PRECISION OF DELAY 0 DELAY RATE 

m 

WT 

!(!!« 



PI COVARIANCE MA'rUIX OF OBSERVABLES 

m 

wr 

D>K 





wr 


WRITE 

s 

SIC A PRIORI VARIANCE OF UNIT WEIGHT 

m 

WT 

m 



PU,P WEIGHT MATRIX OF ODSERVAULES 

Hole 

wr 

.+* 



H1G1.S1G2 

Ht* 

WT 





>11111 

WT 

«!!(; 

ouTi'irr 

! 

P.SIG 


wr 

WT 

4* 

tf* 

OPTIONS 

: 

DELAYS ONLY 

Sfi* 

!|(!K 

wr 

wr 

1H:H 



DELAY RA’IES ONLY 


** 



DELAY AND DEI.AY RA-TES 


WT 





>n* 

wr 



GALLS SSP HUUilOUTlNES LOC,DSINV 


wr 






WT 


IMPLICIT REAL*8( A-If . 0-Z) 

INTEGEU TWOa' l 

DIMENSION P(0. PRd), PHI) 

COMMOM /FLG/ F 1 , I0UM2, 1DUM8, IDUM4/DS/TWU, JU8T 

INPUT WIGin’ING INFORMATION 

SIOl IS THE PRECISION OF TIME DELAY (MEIUHS) 

SIG2 IS THE PHEGISION OF TIME DELAY HATE ( METERS/ HOIHl) 

GALL ERASE 
CALL HOME 
GALL ANMODE 
TOITE («i,2) 

FORMA'l’ (5X, 'INPUT WEIGHTING INFORMATION’//’ SIGl ! PIRICISION OF TI 
INF. DELAY IN MirPEHS’/’ SIG2 : PRF.GISION OF TIME DELAY RATE',’ IN M 
2ETERS/H0UR’ ’ E.G. 0.03 0. 100( CORRESPONDS IX) 0. I NS,0. I PS/S)’) 
READ (3,*) SIG1.S1G2 

KS=Tfm 

IF DFl.AV RATE INCLUDED DODDLE DHffiNSIONS OF WEIGHT MATRIX 
IF(Fl.n0.2) KS«KSM.a 

GOrU’UTE OF ELEMEN'IS IN UPPER TRIANGULAR WEIGHT MATRIX 

KR-TNIW:(TNR+1)/a 

WRITE (0,3) 

FORMAT (/’ INPUT GOVARIANGE MATRIX OF OHSERVA'I’IONS’ /’ ENTER IN UPP 
I EH TRIANGULAR FORM COLUMNWISE - DIAGONAL ELEMEN'IS SCALED TX) UNI’TY’ 
2 ) 

READ (5,!(0 U’K I) , 1=1, KR) 

1F(FI.E0.2) GO TO 7 
1F(F1.E«.3) GO IX) 5 

SCALE FOR DELAY NOISE 
DO 4 KG=1,KR 
P(KG)=PHKG):):SIGI**a 
GO TO H 

SCALE FOR DELAY RATE NOISE 
DO « KC=1,KH 
P( KG) = P 1 < KG) *S1G2**2 
GO TO II 

DI'.VI’.LOP GOVARIANGE MATRIX FOR DELAVRDELAV RATF«S 
KR=KS*(KS+l)/2 


WT 

WT 

WT 

WT 

WT 

WT 

wr 

WT 

WT 

WT 

WT 

WT 

WT 

WT 

WT 

OT 

wr 

wr 

WT 

WT 

WT 

WT 

WT 

w 

WT 

wr 

w 

wr 

wr 

wr 

wr 

WT 

wr 

Wl' 

WT 

WT 

w’r 

WT 

WT 

WT 

WT 

WT 

WT 

wr 


10 

ao 
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70 
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DO 
100 
no 
120 
I no 
140 
ISO 
160 
170 
I HO 
190 
260 
210 
220 
230 
240 
230 
260 
270 
2U0 
290 
300 
310 
320 
330 
340 
3S0 
360 
370 
3B0 
390 
400 
410 
420 
430 
440 
4S0 
460 
470 
4H0 
490 
500 
RIO 
820 
530 
840 
550 
560 
570 
5U0 
890 
600 
610 
620 
630 
640 
680 
660 
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na ft ftft ftft 


9 

10 


1)0 u i>i.Kn 

p( D-o.no 

COMTIWUE 
1 1«0 

DO 10 I-I.TKB 
n«ii+i 
DO 9 .f- 1 , 1 1 
CA1.L LOC <J, 

KA«2*J-l 

KD«2*I“I 

KC»2>l!j 

KD=2*1 

CAM. t,On (KA.KR. ID.KH.KB. 1) 
CALL LOG (KC.KD. IV.KH.KH, 1) 
l'( IT) »S 10 1**2 

P( IV)«P1< IT) *8102**2 
CONTI NUli 
CONTINUE 


INVERT COVAIUANOE NATlllX TO CBT WEIGHT MATRIX 
It OAI.I- nSINV (P.K8, 0.0001, lER) 


810 18 im A PlllOni VARIANCE OP UNIT WEICHT 
RIC«1.D0/P( I) 

8CALE WEIUUT MATRIX 
DO 12 KO-I.KH 
12 P(KC)«P(KC)*SIO 


PHIN’I'OUT WEIOIITINO INFORMATION 
IF(TNB.LE.N) GO TO 13 
CALL ERASE 
OAI L HONE 
CALL ANMODE 

13 WHITE (7,14) 

WRITE (0,14) 

14 EOHMiVT ( // 12X, ’ WEIGUTINO OF OU8ERVATION8 ’ ) 

WRITE (fi.ltl) 

WRITE (7,15) 

13 FORMAT (//9X, ’WEIGHT MATRIX - SCALED TO FIRST ELEMENT DNITV*) 

1 1«0 

DO 20 J«1,KS 
11 = 11+1 
DO Hi Is 1,11 

CALL LOC (J, I, IT,T’HD,TND, 1) 

PIK 1)=P( IT) 

16 CONTINUE 

IF(F1,NE.2) WHITE (6,17) ( l’B( K> ,K= 1 , 1 1) 
lEd'I.KO.H) WHITE (6, IB) ( PB( K) . K» 1 , 1 1 ) 

WHITE (7,19) (PD(K) ,K«1, II) 

17 FORMAT (/,9X,6F7.4) 

IB EOHNAT (9X.6E10.7) 

19 FORMAT (9X, 12F7.4) 

20 CONTINUE 

WRITE (6,21) SIC 

21 FORMAT (/9X, ’A PRIORI VARIANCE OF UNIT WEIGffl’ »’,F10.0) 

WRITE (7,22) 8101,8102,110 

22 FORMAT (//.9X, 'TIME DELi Y’ , 5X,F 10. B, 5X, * (METERS) ’/9X, ’DELAY RATE’ 
15X,F10.5,5X, ’(METEHS/UOIR) ’/9X, ’A PRIORI VARIANCE OF UNIT WEIGUT’ 
2,9X,F10.f>) 

G 

WRITE (6,23) 

23 FORMAT (//' PRESS I - THEN RETURN’) 

READ (3,*) MOM 

C 

CALL ERASE 
CALL ROME • 

CALL ANMODE 
WRITE (6,24) 
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WT 

910 

WT 

920 

WT 

930 

WT 

940 

WT 

950 

WT 

960 

WT 

970 

WT 

900 

WT 

990 

WT 

1000 

WT 

1010 

WT 

1020 

WT 

1030 

wr 

1040 

WT 

1030 

WT 

1060 

WT 

1070 

WT 

tOBO 

Wl' 

1090 

Wl' 

1100 

WT 

1 1 10 

WT 

1120 

Wl' 

1 130 

WT 

1 140 

Wl' 

1130 

wr 

1 160 

WT 

1170 

WT 

1100 

WT 

1 190 

WT 

1200 

WT 

1210 

WT 

1220 

wr 

1230 

WT 

1240 

WT 

1250 

WT 

1260 

WT 

1270 

WT 

1200 

WT 

1290 

WT 

1300 

WT 

1310 

W'l' 

1320 

Wl' 

1330 

wr 

1340 
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EH WElGirriMG DATA ENTER 2 
,• > 

READ (5.SK) NON 
lF(M0n.Ea.2> GO 'lO 1 

RETURI4 

END 


24 


FOHHAT (//, ’ IF YOU WISH TO HEEHT 
1 OTIIEHWIS'E ENTER ANY OTHER NUNHEII 


NT lUSO 
WT 1000 

VT loro 

NT lililO 
VT 1390 
VT 1400 
VT 1410 


oona 


SUBROUTINE FILL ( A.P, VTPVl . AL,B> 

m MOraiAL MATRIX HANDLER ** 

«« 

** FILL NORMAL MATRIX IN A SEttUENTIAL MANNER •** 

** IN TRIANGULAR STORAGE >k* 

«« »» 

** INPUT I A GONIllIBUl'ION TO NORMAL MATRIX AT GIVEN 

m EPOCH OF ALL OBSERVING STATIONS ** 

** V WEIGirr MATRIX OF OBSERVABLES 

** KK TOTAL NUMBER OF PARAMETERS ** 

** FI.F2 PIUIGRAM FLAGS ** 

** AL MISCL08URE VECTOR ** 

** TNB TOTAL NUMBER OF BASELINES t* 

** OUTPUT : V CONTRIBUnON TO NORMAL MATRIX *t. 

m U CONTRIBUTION TO U«ATPL VECTOR ** 

** VTPV CONTRIBUTION TO SUM OF RESIDUALS SOUARED ** 

** 

** CAI.LS SUBROUTINE MATPV, LOC(SSP) ** 

)K!»: 

IMPLICIT REAL*0(A-1I.0-Z) 

INTEGER Fl.FB.TND 

DIMENSION A( 1) . W<1), U(l), P( 1) , B(l), VTPV(l), ALK 12) , AL< 1) 
COMMON /FLG/ FI , IDUMSJ.FB, IDUM4/H'nUI/W/UTR>IXU/BS/'TND,KK 

KS»TNB 

IF(FI.E0.2) KS^KSilcB 
CALCULATE ATP 

CALL MATPV ( A.P.B.KK.KB.KB.O. 1) 

I1»0 
IND» 1 

DO 5 Jsl.KK 
II»II+1 

IFIFO.EQ. 1) GO TO 2 
DO 1 JK«I.KS 

CALL LOC ( J.JK, lY.KK.KB.O) 

ADD CONTRIBUriON TO U»ATPL VECTOR 
U( J)=U<J)-B( IY)*AL(JK) 

CONTINUE 
IN) 4 I»l, II 
DO 3 K»1,KS 
KZ® ( K— 1 ) !(*KK 

ADD CONTRIBUTION TO NORMAL MATOIX 
W( IND)®W( INI))+B( I+KZ)*A( J+KZ) 

CONTINUE 
IND® IND+1 
CONTINUE 
CONTINUE 

IF(F3.Ett. 1) GO TO 6 

ADD CONTRIDUl'ION TO SUM OF RESIDUALS SttUARED 
VTPV=LTPL+XTU 

CALCULAIT LTPL IN THIS ROUTINE 
CALL MATPV ( AL, P, ALl , 1 , KS.KS.O. I ) 

CALL MATPV ( ALl . AI.. VITV, 1 , HS, 1 ,0,0) 

VTPV1®VTPV1+VTPV< 1) 

RETURN 

END 



FL 

10 

:»* 

FL 

20 

m* 

FL 

30 

** 

FL 

40 


FL 

00 


FL 

60 


FL 

70 


FI 

B0 


FL 

90 


FL 

100 


FL 

110 


FL 

120 


FL 

130 

** 

FL 

140 


FL 

130 


FL 

160 


FL 

170 

** 

FL 

180 


FL 

190 


FL 

206 


FL 

210 


FL 

220 


FL 

230 


FL 

240 


FL 

250 


FL 

260 


FL 

270 


FL 

200 


FL 

290 


FL 

300 


FL 

310 


FL 

320 


FL 

330 


FL 

340 


FL 

350 


FL 

360 


FL 

370 


FL 

380 


FL 

390 


FL 

400 


FL 

410 


FL 

420 


FL 

430 


FL 

440 


FL 

450 


FL 

460 


FL 

470 


FL 

480 


FL 

490 


FL 

500 


FL 

510 


FI, 

520 


FL 

530 


FL 

340 


FL 

550 


FL 

560 


FL 

570 


FL 

580 


FL 

590 


FL 

600 


FL 

610 


FL 

620 


FL 

630 
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OOOftO O ft ftftftft ft ftft 


l.JK, IJ.NB) 


** INPUT 

)|c« 

** 

#!)! 

** 

** 

** 

** 

** 

** 

OUTPUT 

** 

** 

*'* OPTIONS 


OUTPUT 


CALCULATES PARTIAL DERIVATIVES 

OHG EARTH ROTATION RATE 

C SPEED OF LIGHT 

ERAD MEAN EARTH RADIUS 

CONV CONVERTS UNIVERSAL TO SIDEREAL TINE 

DK.DY.DZ STATION COORDINATE DIFFERENCES 

RA,D ttUASAR COORDINATES 

IP OUASAR NUMBER 

DS.FRNG OBSERVED DELAY Q DELAY RATE 

JK ONE BASELINE CASE INDEX 

IJ SlEP NUMBER 

NB BASELINE NUMBER 

F 1 , F2 , . . . PROGRAM FLAGS 

ni GST AT INITIAL EPOCH 

TK EPOCH OF OBSERVATION 

DSe.FRNCe THEORETICAL PARTIALS 

AL MI8CL0SURE VECTOR 

PART.G DELAY 8 DELAY RAIE PARTIALS 

DELAYS ONLY 
DELAY BATES ONLY 
DELAY AND DELAY RATES 


IMPLICIT HEAL*8(A-H,0-Z) 

INTEGER F1,F2,F3,F4 

DIMENSION F( 1) , G(l), X(l). Y(l>, Z(l), RA(l), D(l), ALU) 
COMMON /FLG/ F 1 , F8 , F3 , F4/PDR1 /F/'PDR2/'C/'S0U 1 /RA/S0U2/D 

COMPUTE TRIGONOMETRIC MEMBERS OF PARTIAL DERIVATIVES 
CD»DCOS(D< IP)) 

SD=DSIN( D( IP)) 

Yl = OMG*TK+TH 
Y2=Y1-RA( IP) 

CKP»DCOS( Y2) 

SI0*=DSIN(Y2) 

IF(Fl.Ett.a) GO TO 6 

COMPUTE PARTIAL DERIVATIVES 

TAU PARTIAL DERIVATIVE 
F( l)=-CDi)«CKP 

EPSILON PARTIAL DERIVATIVE 
F(2)=CD*SKP 

SIGMA PARTIAL DERIVATIVE 
F(3)=-SD 

EXP= ( DX*SKP+DY*CKP) *CD 
IF( IJ.Ett. I) GO TO 3 

SKIP EARTH ORIENTATION PARAMETERS FOR FIRST STEP 
TO PROVIDE INITIAL ORIENTATION OF BASELINE 
POLAR MOTION DIFFEIUCNGES PARTIAL DERIVATIVES 
SKIP ETA IF F2=3 
IF(F2.Ett.3) GO TO 1 
KSI COMPONENT 
F( 4) =-DX*SD+DZ*CD:ltCKP 
F(4)=F(4)/ERAD 
SKIP KSI IF F2=2 


AL 

PR 

10 


PR 

20 


PR 

30 


PR 

40 


PR 

80 

** 

PR 

60 


PR 

70 


PR 

00 

** 

PR 

90 


PH 

100 

** 

PR 

110 

*)K 

PR 

120 

** 

PR 

130 

** 

PR 

140 


PR 

180 

** 

PR 

160 

** 

PR 

170 


PR 

180 


PR 

190 

** 

PR 

200 


PR 

210 

»» 

PR 

220 


PR 

230 

** 

PR 

240 

** 

PR 

280 


PR 

260 


PR 

270 


PR 

280 


PR 

290 

:»iK 

PR 

300 


PR 

310 


PR 

320 


PR 

330 


PR 

340 


PR 

380 


PR 

360 


PR 

370 


PR 

380 


PR 

390 


PR 

400 


PR 

410 


PR 

420 


PR 

430 


PR 

440 


PR 

430 


PR 

460 


PR 

470 


PR 

480 


PR 

490 


PR 

500 


PR 

510 


PR 

520 


PR 

530 


PR 

540 


PR 

550 


PR 

560 


PR 

370 


PR 

380 


PR 

590 


PR 

600 


PR 

610 


PR 

620 


PR 

630 


PR 

640 


PR 

630 


PR 

660 
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no nn a aa ft ftft ft ft ft ftftft ft ft ft ftft ftft ft ft ftft ft .^oo 


C 

1 


3 


* 

5 


6 

7 


B 


iF(ra.Ea.2> go to a 

PR 

670 

ETA COMPONENT 

PR 

600 

F( 5 ) > DV)«:SD4^ OZ»CD«SKP 

PR 

690 

F(5-JK)«F(B)/E!IA» 

PR 

700 

UTl-int: PAIITIAL DERIVATIVE 

PR 

710 

F(6)«CONV*EXI* 

PH 

720 

CHANGE UNITS TO HILU8ECONDS 

PR 

730 

F(6-JK)“l».D0»F(6)/'( ieOO.D«»RO) 

PR 

740 


PR 

780 

DECLINATIONS 

PR 

760 

F( 7) » ( DX*CKP-DY#SKP) »SD-DZ#CD 

PR 

770 

F(7“JK)=F(7)/R0 

PR 

780 

RIGHT ASCENSION DIFFERENCES 

PR 

790 

F(0)«-EXP 

PH 

800 

F(0-JK)>«F(0)/R0 

PR 

810 


PR 

820 

SKIP CLOCK PARAMETERS IF F4»2 

PR 

830 

IF(F4.Ea.2) GO TO 4 

PH 

840 

CLOCK OFFSET PARTIAL DERIVATIVE 

PR 

880 

F(9-JK)»C 

PR 

060 

CLOCK RATE PARTIAL DERIVATIVE 

PH 

870 

F( I0-JK)»C*TK 

PR 

880 


PR 

890 

SKIP DELAY RATES IF FI^l 
IFIFl.EO. 1) GO TO 9 

PR 

900 

PR 

910 


PH 

920 

DELAY RA'IE PARTIALS 

PR 

900 

OMG«OMG*CD 

PR 

940 

TAU 

PR 

950 

G( l)»OMC*SKI* 

PR 

960 

EPSILON 

PH 

970 

G(2)»OMC*CKl' 

PH 

900 

SIGMA 

PR 

990 

G(3)*O.D0 

PH 

1000 

EXT= ( DX*CKP-DY#SKP) sf^OMC 

PR 

1010 

IF( IJ.EO. 1) GO TO 8 

PH 

1020 

IF(F2,Ett.3) GO TO 6 

PR 

1000 


PR 

1040 

POLAR MOTION DIFFERENCES 

PH 

1030 

KSI COMPONENT 

PR 

1060 

G(4)=-0MC*DZ^!SKP 

PR 

1070 

G(4)=G(4)/EHAD 

PH 

1080 

IF(F2.Ea.2» GO TO 7 

PH 

1090 

ETA COMPONENT 

PH 

1100 

G(n)=OMG*DZ*CKP 

PR 

1110 

G( 5-JK)»C(5)/ERAD 

PH 

1120 

irn-UTC DIFFERENCES 

PR 

1 100 

C( 6 ) a CONVKEXT-I'OMG 

PH 

1 140 

C1L\NGE UNITS TO MILLISECONDS 

PH 

1 150 

G(6-JK)=G(6)/(3600.D0*100O.DO) 

PR 

1 160 


PH 

1 170 

DECLINATIONS 

PH 

1180 

G( 7) =-OMG:t:SD;»:( DX)(cSKP+DY»CKP) 

PR 

1190 

G( 7-JK)»G(7)/RO 

PR 

1200 

RIGin' ASCENSION DIFFERENCES 

PR 

1210 

G( 8) a -EXT 

PH 

1220 

C(Q-JK)«=C(8)/R0 

PR 

1200 

IF(F4.E0.2) GO TO 9 

PR 

1240 


PR 

1250 

CLOCK OFFSET DIFFERENCES 

PH 

1260 

G( 9-JK)=0.DO 

PR 

1270 

CLOCK RATE DIFFERENCES 

PH 

1280 

G( 10-JK)=C 

PH 

1290 


PH 

1000 

Ftf/rURN IF COVARIANCE ANALTOIS ONLY 

PH 

1010 

lF(F;).Ett. 1) GO TO 12 

PH 

1020 


PH 

1000 

COMPUTE APPUOXIMA'IE VALUE OF OBSERVATION 

PR 

1040 
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na 



NC>KB 

PR 

1380 


iFiri.EO. i.0R.ri.Ea.2> go to i« 

PR 

1360 


NC»KC“i 

PR 

1370 


GO TO II 

PR 

1330 

!• 

DiO- >• UX»GD«CKP+OY»GD»SKP>DZ»8D 

PR 

PR 

PH 

1390 

1400 

1410 


CALUULATB HI6CL08UHE8 


IF<F1.EQ.2> NC«2«NB-1 

PR 

1420 


AMItOoDOO-DS 

PR 

1430 


IFd'l.EQ. 1) GO TO 12 

PR 

1440 

11 

FRNG0«ONC»( DJ^OKP'f DV«GKP> 

PH 

1480 


AL( KG+ 1 ) »FHNG0-FRWC 

PR 

1460 


RETURN 

PR 

1470 

12 

PR 

1400 


END 

PR 

1490 
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aa a ct n na ana 


4 


C 

c 


9 


8UBH0UT1NE ANATR (RB, IP, iN.STEP.IfSTBra, JK, IPFIX) 


** 

** 

Hi)|t 

** 

** 

** 

m 

m 

m 

** 

*)K 


DESION HATIilX HANDLER 


INPUT 


: NB 
TWB 
IP 
KX 
IM 

STEP 

KSTEPS 

JK 

IPFIX 

Fl,ra,P4 

F.G 


NUMBER OP BASELINE 
1WAL NUMBER OF BASELINES 
OUASAR NUMBER 
T)TAL NUMBER OF PARAMETERS 
1\)TAL NUMBER OF OUASARS 
EARTH ORIENTATION STEP NUMBER 
1V>TAL NUMBER OF STEPS 
ONE BASELINE CASE INDEX 
ttUASAR OF FIXED RIGHT ASCENSION 
PIIOORAM FLAGS 

PARTIALS OF DELAY 8 DELAY RATE 


OUTPUT s A 


OBSERVATION CONTRIBUTION TO A - MATRIX 


m OPTIONS ! Dl'XAVS ONLY 


!4t« 


DELAY AND DELAY RATES 


IMPLICIT RF,AL*B( A-H.O-Z) 

INTEGER THB,FI.F2,F4,8’nSP 
DIflENSiON F( 1) . G( I) . A( I) 

COMMON /FLC/ F I , F2, IDUMO. F4/PDRI/F/PDR2/G/AT1UL'A/B8/TNB,ICK 
DELAYS 

TAU EPSILON SIGMA 
JI»8*(NB-1) + I 
DO I IJ»l,a 
A( JI)«=F( IJ) 

J1*JI+I 

CONTINUE 

J2»3st!TNB 

SKIP EAR'ni ORIENTATION PARAMEI’EBS FOR FIRST STEP 
•IX) PROVIDE lU:! KRENCE ORIENTATION 
IFIS'IUP.Ea, I) GO TO 3 
POLAR MOTION COMPONENT DIFFERENCES 
J3«J2+STEP“1 
A( J3)«F(41 

IF(F2.E0.2.0R.F2.Ett.3) GO TO 2 
A< J3+NSTEPS-1)»F(5) 

J4» J2+( 2-.IKJ *KS’l'EPS+S'I'EP-( 3-JK) 

UTl DIFFERENCE 
A< J4)»F(6-JK) 

DECLINATIONS 

J0» J2+ ( 3-JK) *( NSTEPS- 1 ) + IP 
A(J5)»F(7-JK) 

IF REFER1-;NCE OUASAR do not FILL DESIGN MATRIX 
ro PROVIDE RIGHT ASCENSION ORIGIN 
1F( IP. Ett. IPFIX) GO TO 4 
J<»= J5+ IM 

IF( IP. GT. IPFIX) J6«J6-I 
RIGirr ASCENSION DIFFERENCES 
A(J6)*F(8-.IK) 

IF(F4.E0.2) GO IX) 5 

J7- J2+( 3-JK) *( NSTEPS- 1 ) +2* IM+2*( NB- 1) 

J0=J7+1 

CLOCK OFFSET DIFFElUiNCES 
A( J7)=F(*)-JK) 

CLOCK MTE DIFFERENCES 
A( JB)=F( 10- JK) 

IFd'l .Ett, i) GO XX) 10 



AM 

IS 

I*ltt 

AN 

20 


AM 

80 

** 

AM 

40 


AM 

80 

** 

AH 

60 


AN 

70 


AM 

00 

** 

AH 

90 

** 

AH 

100 


AH 

110 

** 

AM 

120 


AN 

130 

** 

AM 

140 

** 

AM 

150 

m* 

AM 

160 

** 

AM 

170 

** 

AM 

100 

** 

AM 

190 

** 

AM 

200 

** 

AN 

210 

** 

AM 

„20 

:** 

AM 

230 


AM 

240 


AM 

250 


AM 

260 


AM 

270 


AM 

200 


AM 

290 


AM 

300 


AM 

310 


AM 

320 


AM 

330 


AN 

840 


AM 

350 


AM 

360 


AM 

370 


AM 

300 


AM 

390 


AM 

400 


AM 

410 


AM 

420 


AM 

430 


AM 

440 


AM 

450 


AM 

460 


AM 

470 


AM 

400 


AM 

490 


AN 

500 


AH 

510 


AM 

520 


AM 

530 


AM 

540 


AM 

550 


AN 

560 


AM 

570 


AH 

500 


AH 

590 


AM 

600 


AN 

610 


AN 

620 


AM 

630 


AM 

640 


AM 

650 


AN 

660 
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DELAY RATES 

DO « IJ*1,8 
A<JHKK>>C( IJ> 

JI>JI+I 

COMTINUE 

ir<STEP.EO. 1) €» TO 8 
A<J8't|CK)-G(4> 

IF(F2.E0.2.0R.r2.Ea.8> GO TO 7 
A( J3-»^N81i;PS- 14KK) -Oi 8) 
A(J4+KIO>G(6-JIO 
A(J8‘«-iaO«G(7-JK) 

1F( IP. EG. IPFIX) GO TO 9 
A<J6+nO>G< 8-JllO 
lF(F4.Ea.2> GO TO 10 
A(J74^KK>-G(9-JK) 

A<JB+KK><G( 10- JK) 


RETinUf 

END 


• 4 * 1*4 <4 « 4 *« •««)•«»»» 
•4»5So»—S«e«4»a*ei4— •«s*4 


I 



c 

SUBROUTINE rAfTTR (NB. IP, IM.S'HO'.RSTErfl, JK. imX) 

FM 

FM 

10 

20 

c 


*its 

FM 

30 

c 

** DESIGN MATRIX UANDLER FOR DELAY RATE PARTIALS 

»» 

FM 

40 

c 


m* 

FM 

SO 

c 

** SEE SISBROUTINL AHATR FOR DETAILS 

** 

FM 

<iO 

G 


»)fc 

FM 

70 

C 


FM 

UO 

C 



FM 

90 


IMPLIGIT REAL*0(A-II,G"Z) 


FM 

100 


INTKCiER TNB.F2,F4,8TEP 


FM 

1 10 


DIMENSION 0( 1) . A( I) 


FM 

120 


COMMON /F LG" I DUM 1 , F2 , 1 DUMB , F4/PDR 1 /DUM< 1 ) /PDR2/C^ATRX/A/B8/T»B 

.KK 

FM 

IBO 

C 



FM 

140 

c 

SIGMA, CLOCK OFFSET PARAMETEBS NOT ESTIMABLE 


FM 

ISO 

c 



FM 

160 

c 

TAU EPSILON 


FM 

170 


JI«2*(ND-I)+1 


FM 

lUO 


DO 1 IJ«1,2 


FM 

190 


A(Jn«C( IJ) 


FM 

200 


JI«J1+| 


FM 

210 

1 

CONTINUE 


FM 

220 


J2«2*TND 


FM 

230 

c 

EARTH ORIENTATION PARAMETERS 


FM 

240 


IFiSTEP.Ett. 1) GO TO 8 


FM 

2S0 


J3«J2+STEP~J 


FM 

260 

c 

KBI COMI'ONENT 


FM 

270 


A< J3)=G(4> 


FM 

2110 


IF(F2.Ea,2.0U.Fa.E0.3) CO TO 2 


FM 

290 

c 

ETA COMPONENT 


FM 

300 


A( JO+KSTOPS- I) «G( 5) 


FM 

310 

2 

J4“ J2+ ( 2- JK) *NSTEPS+STEP~( 3-JK) 


FM 

02© 

G 

UTC-Gri 


FM 

330 


A( J4)«G(6-JK) 


FM 

34© 

3 

IF( IP.EQ. IPFIX) CO TO 4 


FM 

33« 


J5= J2+( a- JK) *( N3TEP8-- 1) + IP 


FM 

36© 


IF( IP. GT, IPFIX) J5»J5-1 


FM 

37© 

C 

DECLINATION DIFFERENCES 


FM 

3U© 


A< JS)=G< 7-JK) 


FM 

39© 


J6=.l5+IM-l 


FM 

40© 

C 

RIGHT ASCENSION DIFFERENCES 


FM 

41© 


A< J6>=G(0-JK) 


FM 

42© 

4 

IF(F4.E(1.2) GO TO 3 


FM 

43© 

C 

CLOCK RATE DIFFERENCFJS 


FM 

44© 


J7=J2+(8“JK)*(NSWS- l)+2t( IM-D+NB 


FM 

43© 


A( J7)»G( 10-JK) 


FM 

46© 

C 



FM 

470 

5 

RETLItN 


FM 

4(1© 


END 


FM 

49© 


125 



0 

c 

G 

C 

C 

G 

G 

C 

C 

C 

C 

G 

G 

C 

C 

C 

C 

C 

C 

G 

G 

G 

G 

G 

G 

G 

G 

C 

G 

G 


SUBHOirrilfl iOLVB C«IGt,GOI«,ID{, IH.MIVni.DIIT.B.IHI.IM.VirVI. lOOUIIT BL tf 

n 91 99 


G 

G 

G 


G 

G 


4 

5 


C 

G 

G 


mm 




«« SL 

4S 



LEAST SOUARES ESTIMATION RANDLEM 

m» 

SL 

8t 





mm 

SL 

6« 


INPUT 

t XG,YG,ZG 

BASELINE GONPONENTB 

mm 

SL 

7S 



N 

COMPLETE NORMAL MATRIX 

mm 

SL 

00 

«« 


U 

COMPLETE ATPL VECTOR 

mm 

SL 

90 

mm 


VTPVi 

LTPL CONTRIBUTION TO VTPV 

mm 

SL 

too 

mm 


S1(XI 

A PRIORI VARIANCE OP UNIT WBIOn' 

mm 

SL* 

110 

mm 


DiST 

BASELINE DISTANCES 

mm 

SL 

120 

mm 


KK 

TOTAL NUMBER OF PARAMETERS 

mm 

SL 

180 

mm 


IN 

total number of quasars 

ma 

SL 

140 

mm 


TNB 

TOTAL NUMBER OP BASELINES 

mm 

SL 

180 

mm 


NSTEP8 

TOTAL NUMBER OF STEPS 

mm 

SL 

160 

mm 


FI. . . 

PROGRAM FLAGS 

mm 

SL 

170 

mm 


ICOUNT 

TOT'AL NUMBER OF OBSERVATIONS 

mm 

SL 

180 

mm 




mm 

81. 

190 

mm 

WRITE 

i N 

NORMAL a COVARIANCE MATRIX 

mm 

SL 

800 

mm 


DM 

BASELINE VARIANCE 8 COVARIANGES 

mm 

SL 

210 

mm 


CORR 

CORRELATION MATRIX OF PARAMETBS 

mm 

SL 

220 

mm 


VTTV 

SUM OF RESIDUALS SQUARED 

mm 

SL 

230 

mm 


SIG2H 

A POSTERIORI VARIANCE OF UNIT WBIGBT 

mm 

SL 

240 

mm 


XX 

CORRECTION TO APPROXIMATE PARAMETERS 

mm 

SL 

280 

mm 


EIG 

EIGENVALUES OF COVARIANCE NATIIIX 

mm 

SL 

260 

mm 




mm 

SL 

270 

mm 


GALLS 88P 

ROUTINES : I^INV.DBIGEN.LOC 

mm 

SL 

280 

mm 


CALLS SUBROUTINES STDLST 

mm 

SL 

290 

mm 




mm 

SL 

800 




IMPLICIT Iffi,AL*0(A-H,L-Z) 
IN'nSGER TNB,Fl,r2.F8,F4,NSTBPS 
DIMENSION }{G( 1) . YG<1), ZC(I), 
IGORRd), lUl). XX(1> 


DIST(l), EM(1>, DN(1>, BID, N(l>, 


SL 

SL 

SL 

SL 

SL 

SL 


COMMON /FLG/ F 1 , F2 , F3 , F4/CHD4/XC/CRDB/'YC/CBD6/'ZC/NTH)^N/0TIIX/O/BS/' SL 


nifB^KK 

PRINT NUMBER OF OBSERVATIONS 
DOUBLE OBSERVATIONS IF DELAY RATE INCLUDED 
IF( F 1 . EO. 2) ICOUNT* ICOUNTVi2 
WRITE (7 1) ICOUNT 

FORMAT (//I2X, 'THE NUMBER OF OBSERVATIONS «M4) 

PRINTOUT NORMAL MATRIX 
WRITE (7,2) 

FORMAT ( IHl, I2X, ’NORMAL MATRIX’) 

II»0 

DO 5 J*1,KK 
11*114^1 
DO 3 let, II 

GALL LOG (J, I, IT.KK.KK, 1) 

CORR( I)»N( IT) 

CONTINUE 

WRITE (7,4) J,(COIUl(K),KaI,II) 

FORMAT (/,8X, 12,’. ’ , (TIO, I2F9.2) ) 

CONTINUE 

INVERT NORMAL MATRIX 

ONLY UPPER TRIANGULAR PART IS REEDED 

CALL DSINV (N.KK.e.OOei, lER) 

PRINTOUT VARIANCE-COVARIANCE HA'IRIX (UNSCALED) 

WRITR (7,6) 

FORMAT ( lUI , I2X, ’VARIANCE-COVARIANCE MATRIX (URSCALED) 
II-O 


SL 

SL 

SL 

SL 

SL 

SL 

SL 

SL 

SL 

SL 

SL 

SL 

S,. 

SL 

SL 

SL 

SL 

SL 

SL 

SL 

SL 

SL 

SL 

SL 

SL 

SL 

SL 

SL 

SL 


310 

320 

330 

340 

380 

360 

870 

300 

390 

400 

410 

420 

430 

440 

480 

460 

470 

400 

490 

800 

810 

820 

830 

840 

380 

860 

870 

800 

890 

600 

610 

620 

630 

640 

680 

660 


126 


fto no on 


a 

9 


!• 


C 


It 

12 


18 

14 

18 


16 

17 

la 


19 


C 

20 


190 9 J«M0C 

8L 

670 


SL 

600 

DO 7 l»l,li 

8L 

690 

CALL LOC (J.l, (T.KK.KK, 1) 
COlUU l)»m IT) 

SL 

700 

SL 

710 

coirriNDC 

SL 

720 

MUTE (7,0) J,(COitR<K),iO>l,II) 
POHMAT (/.8X. 12,*. ' , (TIO. I2P9.8) ) 

SL 

780 

SL 

740 

COMTIIlUB 

SL 

780 


SL 

760 

COMPOTE AND LIST STANDARD DEVIATIONS (A PRIORI) 

SL 

77© 

CALL STDLST (SIC^.NSTEPS, IM.XX, 1) 

8L 

700 

SL 

790 

IF DELAV RAIT^- DISTANCE NOT ESTIMABLE 

SL 

000 

IF(FI.Ett.a) 00 TO 01 

St 

010 


SL 

020 

C0NPU1*E DISTANCE COVARIANCE NA11UX 

SL 

030 

J«U»TNE*'n»B 

SL 

040 

DO 10 l«l,J 

SL 

080 

EN( 1)»0.D0 

SL 

860 

U( |)»0.D0 

SL 

070 

JL»0*TNB 

SL 

800 

DO 11 K'l.TTID 

SL 

090 

JK»8*<K-l)*(TNA+l) + l 

SL 

900 

coMPim; Rmum propagation partials 

SL 

910 

B( JK)«XC( K) /DI8T(K) 

SL 

920 

B< JK+1)«VC( Kl/OISTdO 

SL 

930 

U( JK+2) «7C(K)/DIST( K) 

SL 

940 

CONTINUE 

SL 

980 

IB2»TRD*THB 

SL 

960 

DO 12 I>l. 1D2 

SL 

970 

0M( 1)00. DO 

SL 

900 

DO 18 1*1, TNB 

SL 

990 

I8»< I-1)*JL 

SL 

1000 

DO 14 K»1,JL 

SL 

1010 

12-I8+K 

SL 

1020 

DO 13 J»1,JL 

SL 

1030 

CALL LOC (J.K, in,KK,KK, 1) 

SL 

1040 

I4« I8+J 

SL 

1030 

EMC i2)»BN( I2)+B( I4)«N( IR)^81G2 

SL 

1060 

CONTINUE 

SL 

1070 

CONTINUE 

SL 

1000 

CONTINUE 

SL 

1090 

K3»0 

SL 

1100 

DO 18 I«1,TNB 

SL 

1110 

I3»( I-1)*JL 

SL 

1120 

DO 17 K«1,TNB 

SL 

1 130 

K2»{K-I)*JL 

SL 

1 140 

K3«K3+I 

SL 

1180 

DO 16 Jxl.JL 

SL 

1 160 

Jl» 18+J 

SL 

1170 

J2“Kn+J 

SL 

1100 

DN( K3) » DM( K3) +EM( J 1 ) *D( J2) 

SL 

1 190 

CONTINUE 

SL 

1200 

CONTINUE 

SL 

1210 

CONTINUE 

SL 

1220 

CALL ERASE 

SL 

1230 

CALL HOME 
CALL ANMODE 

SL 

1240 

SL 

1250 

MUTE (7,19) 

SL 

1260 

FORMAT ( I HI.///, 12X,’ DISTANCE COVARIANCE MATRIX*/) 

SL 

1270 

DO 22 I«1,TNB 

SL 

1200 

DO 20 Jsl.TNB 

SL 

1290 

CALL LOC ( I.J, IR,TNB,'rNB,0> 

SL 

1300 

C1I,\NCE UNITS TO CENTIMETORS SQUARED 

SL 

1310 

CORR(.r)»DM( IR)iC 10000. DO 
CONTINUE 

SL 

1320 

SL 

1330 

MUTE (7.21) (COlUWK) ,K«1,TNB) 

SL 

1340 


127 


ai 

aa 


ao 


24 

ao 


24 


DeVIATWHM 


27 

aa 

29 

s 

so 


c 

c 


81 

32 


S3 

34 

38 


C 

C 


C 

C 


36 


C 

C 


roRHAT ( lox.erio.s) 

COKTINOE 

PHtHT BAHELINE STAJIDARO 
WRITE ( 6 »a 8 > 

WRITE < 7.231 

roimT < // lax. * baseline stanbarb deviation (om * > 

DO as l»l,TNB 

CALL LOC ( I , I , IR.TND.'ITfBLO) 

CHANCE UNITS TO CENTinETEfIS 
BSTD«DSaRT< DN< IR) >«IOO.DO 
imiTE ( 6 . 34 ) I.BSTD 
WRITE < 7 , 24 ) l.BSTO 
FORMAT </, I 2 X, 12 , * . * , IX,F 9 , 8 ) 

CONTINUE 

COMPUTE BASELINE CORRELATION HA'IRIX 
WRITE < 6 , 26 ) 

WRITE (T 

FORMAT <//'/',l 2 X, ’BASELINE CORREUTIOR MATRIX’ I 
Il«l 

DO 29 I«l, 11 fB 
DO 27 J» 1 . 1 1 

CALL LOG < I. 1 , IR.TNB.TNB.O) 

CALL LOC (J.J, IS.TRB.'IRB.O) 

CALL LOC < I,J, IT,TNB,TRB, 0 > 

COmU DN< IT) /DSORTi DN< IR)»DM( IS) > 

CONTINUE 

WHITE < 6 , 28 > 1 ,<C 0 RH(K),K> 1 ,II> 

HRI'lE < 7 , 20 ) I,(CORR<K).K«t,!n 
FORMAT </, I 2 X, 12 , ’, 6 F 8 . 2 ) 

1 I«II +1 

CONTINUE 

WRITE < 6 , 30 ) 

FORMAT <//', 6 X, ’PRESS 1 THEN RETURN’) 

READ ( 0 ,«) 

COMPUTE PARAMETER CORRELATION MATRIX 
WHITE ( 7 , 32 ) 

FOHMAT < lUl, I 0 X, ’PARAMETER CORRELATION MATRIX’//) 

1 I «1 

DO 3 S 1 - 1 . KK 
m 33 J*l, II 

CALL LOG (I, I, m.KK.KX, S) 

CALL LOC (J,J, IS.KK.KK, 1 ) 

CALL LOC ( I,J, IT.KK.KK, 1 ) 

COlUW J) »N( IT) /DSORT( N( IR) »N< IS) ) 

CONTINUE 

WRITE ( 7 , 34 ) I,(COHR(K) ,K- 1 , ID 
FORMAT (/, IX, 12 , ’ . ’ ,(T 6 , 28 F 5 . 2 )) 

II-Il+l 

CONTINUE 

SKIP SOLUTION IF INTERESTED ONLY IN COVARIANCE ANALYSIS 
IFCFO.Ett. I) GO W 41 

CALL ERASE 
CALL HOMF, 

CALL ANMODE 

COMPUTE CORRECTIONS TO APPROXIMATE PARAMETERS 
DO »6 I»I,KK 
DO 36 ,r-i,KK 

CALL LOC < I,.I,m,KK,KK, t> 

XX( I)»XH( I)+N( 1 R)*U(J) 

CONTINUE 

CALCUibATE VTPV 
VrPV 2 »O.D 0 


SL 1380 
8L 1360 
SL 1370 
SI. 1300 
SL 1390 
SL 1400 
SL I4i0 
SL 1420 
EL 1430 
SL 1440 
SL 14S0 
SL 1460 
SL 1470 
SL 1400 
SL 1490 
SL IS00 
SL 1510 
SL 1020 
SL 1830 
SL 1840 
SL 1530 
SL 1860 
SL 1070 
SL 1SH0 
SL 1890 
SL 1600 
SL 1610 
SL 1620 
SL 1630 
SL 1640 
SL 1680 
SL 1660 
SL 1670 
SL 1600 
SL 1690 
SL 1700 
SL 1710 
SL 1720 
SL 1730 
SL 1740 
SL 1780 
SL 1760 
SL 1770 
SL 1700 
SL 1790 
SL 1000 
SL 1010 
SL 1020 
SL 1030 
SL 1040 
SL 1080 
CL 1060 
SL 1070 
SL 1000 
SL 1090 
SL 1900 
SL 1910 
SL 1920 
SL 1930 
SL 1940 
SL 1980 
SL 1960 
SL 1970 
SL 1900 
SL 1990 
SL 2000 
SL 2010 
SL 2020 
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DO 8? I>1,KK 
VTPV2«VTPVa»XXf I)*0( I) 

SL 2080 
SL 2040 

37 

COMTIMUE 

SL 2050 


VTPVF-VTPVl+VTPVa 

SL 2060 

C 


SL 2070 

C 

GALGULATE THE A POSTERIORI VARIAHCE <MP OHIT MBIGVr 

SL 2080 


SIC2H»VIPVr/'( ICOUKT-IOC) 

SL 2090 

G 

CALCULATE A POSTERIORI /A PRIORI 

SL 2100 


6IGR«S!G2H/SIG2 

SL 2110 

G 


SL 2120 

G 


SL 2130 

G 

COMPUTE 3 LIST STAHDARD DEVIATIONS - A POS1VRIORI 

SL 2140 


CALL 8TDLST (Slf^H.MSTEPS. 1H.XX,2> 

SL 2100 

G 

GALL ERASE 

SL 2160 
SL 2170 


GALL HOME 

SL 2180 

G 

CALL ANNODE 

SL 2190 
SL 2200 


WRITE (6,38) VITVP 

SL 2210 


WRITE (7,38) vrpvr 

SL 2220 

aa 

FORMAT (//. 12X, ’VTPV -*,D18.8) 

SL 2230 


WRITE (6,39) SIG2H 

SL 2240 


WRITE (7,39) SIG2H 

SL 2250 

39 

FORMAT (/,I2X, 'A POSTERIORI VARIAHCE OP DHIT WEIGVT -*,D18.8) 

SL 2260 


WHITE (6,40) SIGR 

SL 2270 


WHITE (7,40) SIGR 

SL 2280 

40 

FORMAT (/,12X, ’A POSTERIORI/A PRIORI >*,D10.8> 

SL 2290 

G 

WRITE (6,30) 

SL 2300 
SL 2310 


HEAD (8,») MOM 

SL 2320 

G 


SL 2330 

G 

COHIPUTE EIGENVALUES OF GOVARIAHGE MATRIX 

SL 2340 

41 

WRITE (7,42) 

SL 2t»'0 

42 

FORMAT ( IHl, 18X, 'EIGENVALUES V) 

SL 2 oO 


CALL DEIGEN (N,R,KK, 1) 

SL 1^370 


DO 44 I»1,KK 

SL 2380 


CALL LOG (I,I,IND,KiC,KK, 1) 

SL 2390 


EIG»N( IND) 

SL 2400 


WRITE (7,43) EIG 

SL 2410 

43 

FORMAT (15X,D14.3) 

SL 2420 

44 

CONTINUE 

SL 2430 

G 

GALL ERASE 

SL 2440 
SL 2450 


CALL HOME 

SL 2460 


CALL ANMODE 

SL 2470 


RETURN 

SL 2480 


END 

SL 2490 
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fio a ft ftft ftft 


C 

c 

G 

G 

G 

G 

G 

G 

G 

G 

G 

G 

G 

G 

G 

G 

G 

G 

G 

G 

G 


SUBROUTINE 8TDLST (SIG2.N8TEPS, IH.XX, I8TD) SD 


GOMPUTE a OUTPUT E8TINATE0 STANDARD DEVIATIONS 


2 

3 

4 

C 

5 

6 

7 

8 


«i(C 

** INPUT I 

)tc)K 

** 

** WRITE i 8IGHA 
** OPTIONS 


N 

SIG2 

ri... 

NSTEPS 

TNB 

KK 

IN 

XX 

ISTD 


GOVARIANGE MATRIX OF PARAMETERS 

VARIANGE OF UNIT WEIGHT 

PROGRAM FLAGS 

TOTAL NUMBER OF SIICPS 

TOTAL NUMBER OF BASELINES 

TOTAL NUMBER OF PARAMETERS 

1X)TAL NUMBER OF OUASARS 

PARAMETER GORREGTIONS VECTOR 

INDEX - A PRIORI OR A POSl'BRIORI STD’S 

STANDARD DEVIATIONS OF PARAMETERS 


OUTPirr A PRIORI OR A POSTERIORI STANDARD DEV. 
OUTPUT CORRECTIONS TO APPROXIMATE PARAMETERS 


»He SD 
IK* SD 
«Mc SD 
stent SD 
stene SD 
sleste SD 
stent SD 
nene SD 
nestt SD 
nene SD 
nene SD 
nene SD 
nene SD 
nene SD 
nene SD 
nene SD 
nene SD 
nest: SD 


^nenEnenentnenestests.tnisttntnenentstcnenenentnenenentnenenenentnestenenenenentnenenentntntnenentnenentntnentntneikntntttenknencntntnen: SD 

SD 

IMPLICIT REALneO<A-H,L~Z> 

I N'PEGER NPARM, TNB , F 1 , F2 , F3 , F4 , NSTEPS , TNP 
DIMENSION N( 1) . XX< 1) 

COMIKIN /nTKK/ N/B8/TNB,KK/FLG/F1 ,F2,F8,F4 


STANDARD DEVIATION FUNCTION 
8(A)>DSCUlT(AneSIG2> 

NFaRM is the number OF PARAMETERS COUNTER 
NPARN>e 

IFLOWrO 

IF(S1G2.GT. l.D-10) GOTO I 

SIGNA'O.DO 

IFLOW= 1 

CALL ERASE 
CALL HOME 
GALL ANMODE 


COMPUTE STANDARD DEVIATIONS 
IF( ISTD.E0.2) GO TO 3 
WRITE (6,2) 

WRITE (7 2) 

FORMAT ( i HI, 5X, ’STANDARD DEVIATIONS - A PRIORI’) 

GO IX) 5 
WRITE (6,4) 

WRITE (7,4) 

FORMAT ( IHI,5X, ’STANDARD DEVIATIONS - A POSTERI OR! ’ , / , BX, ’ + PARAHB SD 
ITER CORRECTIONS’) SD 

SD 

IF(Fl.Ea.3) GO TO 7 SD 

WRITE (6,6) SD 

WRITE (7,6) SD 

FORMAT (/,7X, 'TAU EPSILON SIGMA (CM)’) SD 


SD 

SD 

SD 

SD 

SD 

SD 

SD 

SD 

SD 

SD 

SD 

SD 

SD 

SD 

SD 

SD 

SD 

SD 

SD 

SD 

SD 

SD 

SD 

SD 

SD 

SD 

SD 

SD 


IT»TNBne3 
TNP=3 
GO IX) *) 

WRITE (6,8) 

WRITE (7,8) 

FORMA T ( / , 7X, ’ TAU EPS I LON 

IT=TNUnt2 

TNP“2 

ICeO 


(CM) ’) 


SD 

SD 

SD 

SD 

SD 

SI) 

SD 

SD 

SD 


10 

20 

30 

40 

90 

60 

70 

80 

90 

100 

110 

120 

130 

140 

150 

160 

170 

180 

190 

200 

210 

220 

230 

240 

290 

260 

270 

280 

290 

300 

310 

320 

330 

340 

390 

360 

370 

380 

390 

400 

410 

420 

430 

440 

450 

460 

470 

480 

490 

900 

510 

920 

930 

540 

990 

960 

970 

580 

990 

600 

610 

620 

630 

640 

690 

660 


130 



DO 15 I-l.TTIB 

SD 

670 


VWTE (6, 10) I 

8D 

6B0 


WUTE (7, 10) I 

SD 

690 

19 

FORMAT (9X. ’BASELIlfE * M2) 

SD 

700 


DO 14 J»1,TNP 

SD 

710 


IC« 1C+ 1 

SD 

720 


NPARM»NPARM+ 1 

SD 

730 


IF< IFLOW.Eft. 1) GO TO 12 

SD 

740 


CALL LOG ( IG, IG, lO.KK.KK, 1) 

SD 

750 


6IGHAbS(N( 10)) 

SD 

760 

G 

CHANGE UNITO TO CENTIMETERS 

SD 

770 


SIGHA>»SIGHA»100.DO 

SD 

780 


1F( ISTD.E0.2) GO TO 12 

SD 

790 


WRITE (6,11) NP ARM, SIGMA 

SD 

800 


WRITE (7,11) NPARM, SIGMA 

SD 

810 

11 

FORMAT <4X. 13, ' . ’ , 1X,FB.3) 

SD 

820 


GO TO 14 

SD 

830 

12 

XIC-XX( IG)9|cl00.D0 

SD 

840 


WRITE (6,13) NPARM,SIGMA,XIG 

SD 

850 


WRITE (7,13) NPARM, SIGMA.XIC 

SD 

860 

13 

FORMAT (4X, 13, ’ . ’ , 1X,FB.3,3X,F12.4> 

SD 

870 

14 

CONTINUE 

SD 

880 

IS 

CONTINUE 

SD 

890 

C 


SD 

900 


WRITE (6,16) 

SD 

910 


WRITE (7, 16) 

SD 

920 

16 

FORMAT (/,7X, ’POLAR MOTION VARIATIONS (CM)’) 

SD 

930 


IF(F2.E0.2) WRITE (6,17) 

SD 

940 


IF(F2.E0.2) WRITE (7,17) 

SD 

950 

17 

FORMAT (/,9X, ’FIRST COMPONENT ONLY’) 

SD 

960 


IF(F2.E0.3) WRITE (6,18) 

SD 

970 


IF(F2.E0.8) WRITE (7,10) 

SD 

980 

IB 

FORMAT (/,9X, ’SECOND COMPONENT ONLY’) 

SD 

990 


IT=IT+1 

SD 

1000 


ITM= IT+2»( NSTEP8- 1) - 1 

SD 

1010 


IZ=IT+( lTM-IT)/2+l 

SD 

1020 


IF( F2 . EO. 2 . OR. F2 . EO. 3) ITM* lT+NSTEPS-2 

SD 

1030 


DO 22 I=IT, ITM 

SD 

1040 


IF(F2.E0. 1. AND. I.Ett. IT) WRITE (6,19) 

SD 

1050 


IF(F2.Eft. l.AND. I.EO. IT) WRITE (7,19) 

SD 

1060 

19 

FORMAT (9X, ’FIRST COMPONENT’) 

SD 

1070 


IF(F2.Ea. l.AND. I.EO. IZ) WRITE (6,20) 

SD 

1080 


lF(F2.Ea. (.AND. I.EO. IZ) WRIIE (7,20) 

SD 

1090 

20 

FORMAT (9X, ’SECOND COMPONENT’) 

SD 

1100 


NPARM=NPARM+1 

SD 

1110 


IF( IFLOW.Ett. 1) GO TO 21 

SD 

1120 


GALL LOO ( I, I, ia,KK,iCK, 1) 

SD 

1130 


SIGMAsS(N( lO)) 

SD 

1140 

C 

CHANGE UN ns TO CENTIMETERS 

SD 

1150 


S1GMA»S1GMA«1O0.DO 

SD 

1160 


1F( ISTD.Ett.2) GO TO 21 

SD 

1170 


WilTE (6,11) NPARM, SIGMA 

SD 

1180 


WRITE (7,11) NPARM, SIGMA 

SD 

1 190 


GO TO 22 

SD 

1200 

21 

XIC»XX( 1)^100.D0 

SD 

1210 


WRITE (6,13) NPABM,SIGMA,XIC 

SD 

1220 


WHITE (7,13) NPARM, SIGMA,XIC 

SD 

1230 

22 

CONTINUE 

SD 

1240 

C 


SD 

1250 


WRITE (6,23) 

SD 

1260 


WRITE (7,28) 

SD 

1270 

23 

FORMAT (/,7X, ’UTl-UTC VARIATIONS ( 10««2 HICROSEGS) ’ ) 

SD 

1280 


IT=ITMil 

SD 

1290 


ITM= IT+NSTEPS“2 

SD 

1300 


DO 25 1=IT, ITM 

SD 

1310 


NPARM=NPARM+1 

SD 

1320 


1F( IFLOW.Ett. 1) GO TO 24 

SD 

1330 


CALL LOG (1,1, Itt, KK, KK, 1 ) 

SD 

1340 
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SIGHA-SUK HD) 

8D 

1350 

0 

CHANGE UNITS TO 10«»2 NIG!I08EC0ND8 

8D 

1360 


SICHA»SIGHA«ie.De 

8D 

1370 


1F< 1STD.E0.2) GO TO 24 

8D 

1380 


VHITE (6,11) NFARH, SIGMA 

SD 

1390 


WRITE (7.11) NPARM.8IGHA 

SD 

1400 


GO TO 25 

SD 

1410 

24 

xic*xx( DiMe.oe 

SD 

1420 


WHITE (6,13) NPARH,SIGNA,X1G 

SD 

1430 


WHITE (7,13) NPAHM,8IGNA,X.^G 

SD 

1440 

2B 

CONTINUE 

SD 

1450 

G 


SD 

1460 


IF(F(,Ea.3) GO TO 27 

SD 

1470 


WRITE (6,26) 

SD 

1480 


WRITE (7,26) 

SD 

1490 

26 

FORMAT (/,7X, ’DECLINATIONS (NILLIADOSEGS) *) 

SD 

1500 


GO TO 29 

SD 

1510 

27 

WRITE (6,28) 

SD 

1520 


WRITE (7,28) 

SD 

1530 

28 

FORMAT (/,7X, ’DECL. DIFFERENCES (MILLIARC8BC8) * > 

SD 

1540 

29 

IT= ITM+1 

SD 

1550 


ITIM*IT+IM-1 

SD 

1560 


IF(F1.E0.3) ITIM«ITIM-1 

SD 

1570 


DO 31 I=IT, ITIM 

SD 

1580 


NPAHM*NPARM+1 

SD 

1590 


IF( IFLOW.Ett. 1) GO TO 30 

SD 

1600 


GALL LOG (1,1, Itt,KK,KK, 1) 

SD 

1610 


SIGMA»S(N( 10)) 

SD 

1620 

C 

CHANCE UNll'S TO MILLISECONDS 

SD 

1630 


S IGMA=S IGNAiK 1000. D0 

SD 

1640 


IF( 18TD.E0.2) GO TO 30 

SD 

1650 


WRITE (6,11) NPARM, SIGMA 

SD 

1660 


WRITE (7,11) NPAHM, SIGMA 

SD 

1670 


GO TO 31 

SD 

1680 

30 

XIG-XX( I)«10OO.DO 

SD 

1690 


WRITE (6,13) NPARM, SIGMA.XIC 

SD 

1700 


WRITE (7,13) NPAHM, SIGMA.XIC 

SD 

1710 

31 

CONTINUE 

SD 

1720 

C 


SD 

1730 


WRITE (6,32) 
WRITE (7,32) 

SD 

1740 


SD 

1750 

32 

FORMAT (/,7X, ’R.A. DIFFERENCES (NILLIAI1CSEG8) • ) 

SD 

1760 


ITM=ITIM+1 

SD 

1770 


ITN= ITIM+ IM- 1 

SD 

1780 


DO 34 1= ITM, ITN 

SD 

1790 


NPARK=NPARM+1 

SD 

1800 


IF( IFLOW.EO. 1) GO TO 33 

SD 

1810 


GALL LOG ( I, I, lO.KK.KK, 1) 

SD 

1820 


8IGMA»S(N( 10)) 

SD 

1830 

C 

CHANGE UNITS TO MILLISECONDS 

SD 

1840 


S IGMAsS IGMA« 1000 . DO 

SD 

1850 


IF( 1STD.E0.2) GO TO 33 

SD 

I860 


WRITE (6,11) NPARM, SIGMA 

SD 

1870 


WRITE (7,11) NPARM, SIGMA 

SD 

1880 


GO TO 34 

SD 

1890 

33 

XIG«XX( 1) #1000. DO 

SD 

1900 


WRITE (6,13) NPARM, SIGMA.XIC 

SD 

1910 


WRITE (7,13) NPARM, SIGMA.XIC 

SD 

1920 

34 

CONTINUE 

SD 

1930 

C 


SD 

1940 


IF(F4.E0.2) GO TO 45 

SD 

1950 


IF(F1.E0.3) GO TO 39 

SD 

1960 


I F( NPARM.lt. 42) CO TO 35 

SD 

1970 


WRITE (6,46) 

SD 

1980 


READ (9,#) MOM 

SD 

1990 


CALL ERASE 

SD 

2000 


CALL ROME 

SD 

2010 


GALL ANMODE 

SD 

2020 
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G 


SD 

2080 

85 

WRITE (6,86) 

SD 

2040 


WRITE (7,86) 

SD 

2050 

86 

FORMAT (/.rX.'GLOGK OFFSET (IfSEGS)') 

SD 

2060 


ITM«8»TNB+3*( NSTEPS- 1 ) +2» IN 

SD 

2070 


IF(F2.Ea.2.0R.F2.E0.8) 1TN« ITIM M8TEP8-1) 

SD 

2080 


ITN« ITM+2*TNB--2 

SD 

2090 


DO 88 I»ITM, ITN,2 

SD 

2100 


WPARM-NPARM+l 

SD 

2110 


IF( IFLOW.EO. 1) GO TO 87 

SD 

2120 


GALL LOG ( I, I, la.KK.KK, 1) 

SD 

2130 


SIGMA«8(N( 10) ) 

SD 

2140 


IF( 1STD.E0.2) GO TO 87 

SD 

2150 


WRITE (6,11) NPARN.SIGNA 

SD 

2160 


WRITE (7,11) NPARM, SIGMA 

SD 

2170 


GO TO 30 

SD 

2180 

87 

WRITE (6,13) NPARM,SIGMA,XX( I) 

SD 

2190 


WRITE (7,13) MPARM,SIGMA,XX( I) 

SD 

2200 

38 

CONTI HUE 

SD 

2210 

G 


SD 

2220 

89 

WRITE (6,40) 

SD 

2230 


WRITE (7,40) 

SD 

2240 

40 

FORMAT ( / , 7X, ' CLOCK MTE ( P IG0SEC8/HR) * ) 

SD 

2250 


IF(F1.EQ.8) GO TO 41 

SD 

2260 


ITM1“ITM+1 

SD 

2270 


rraisiTN+i 

SD 

2200 


JJ=2 

SD 

2290 


GO TO 42 

SD 

2300 

41 

ITMl* ITN+1 

SD 

2810 


ITN1«ITM1+TNB-1 

SD 

2320 


JJ=1 

SD 

2830 

42 

DO 44 1= iTMi, mn,jj 

SD 

2340 


NPARM- NPARM+1 

SD 

2350 


IF( IFLOW.EO. 1) GO TO 43 

SD 

2360 


CALL LOG (1,1, IO,KK,KK, 1) 

SD 

2370 

G 

CHARGE UNITS TO PIC0SEC8/HR 

SD 

2880 


SIGMA-S(N( IO))«1OO0.D0 

SD 

2890 


IF( ISTD.E0.2) GO TO 48 

SD 

2400 


WRITE (6,11) NPARM, 8 IGMA 

SD 

2410 


WRITE (7,11) NPARM, 8 IGMA 

SD 

2420 


GO TO 44 

SD 

2430 

43 

XCI“XX( I)*1000.D0 

SD 

2440 


WRITE (6,13) NPARM, 8 IGMA, XIC 

SD 

2450 


WRITE (7,13) NPARM, SIGMA, XIC 

SD 

2460 

44 

CONTINUE 

SD 

2470 

G 


SD 

2480 

45 

WRITE (6,46) 

SD 

2490 

46 

FORMAT (//6X, 'PRESS 1 THEN RETURN’) 

SD 

2500 


READ (5,#) MOM 

SD 

2510 

G 


SD 

2520 


RETURN 

SD 

2530 


END 

SD 

2540 
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I 
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SUBROUTINE HAD ( IDEG.NIIV. SEC. AMGLE.PI) RD 

«« »» HD 
** CONVERl>} ANGLE IN DEGREES, MINUTES, SECONDS TO RADIANS «» RD 
** *♦ RD 

RD 

IMPLICIT RE.VL*IKA-U,0-Z) RD 
C*0.D0 RD 
A>SEC/36DO.De RD 
D>HIN/60.De RD 
IF( IDEO.LT.O) G«A+H>IDEG RD 
C'**-C RD 
IKUDEG.GE.O) CA+B+IDEG RD 
ANGLE* C RD 
ANGLE* ANGLE»PI/ 100. DO RD 
lUmiRN RD 
ICND RD 
SUrmoUTINE DEGMS (ANGLE.PI, ideg.hin.seo dg 

** *# DC 
** CONVEHIS ANGLE IN RADIANS TO DEGREES, Ml NinES, SECONDS ** DG 
** ** DG 


IMPLICIT REAL^UI A-II.O-Z) 
ANGLE* ANGLE* 1 80 . UO/P 1 
IDEG* IDINTI ANGLE) 

A«I)FLOAT( I DEG) 

Ct^ ANGLE- A 
(>C*«O.DO 
MIN»IUINT(C) 

B«UrLOAT(MIN) 

I)»C-D 

SEC»I)*60.D0 

D8EC« SEC-60. DO 

DABSEC«0AB8( DSEC) 

1F(D<VD8EC.GT. l.D-9) GO TO 1 

SEC>O.DO 

MIN»M1N+I 

CONTINUE 

1F(M1N.LT.60) GO TO 2 
HIN# MIN-60 
IDEG*IDEG+1 
CONTINUE 

ANGLE* ANGLE*P 1/ 100 . DO 

Ri'/rUHN 

END 


DG 

DG 

DC 

DG 

DG 

DG 

DG 

DG 

DG 

DG 

DG 

DG 

DG 

DG 

DG 

DC 

DG 

DG 

DG 

DG 

DG 

DG 

DG 

DG 


10 

20 

30 

40 

30 

60 

70 

00 

90 

100 

no 

120 

130 

140 

150 

160 

170 

100 

10 

20 

30 

40 

50 

60 

70 

00 

90 

100 

no 

120 

130 

140 

150 

160 

170 

100 

190 

200 

210 

220 

230 

240 

250 

260 


270 

200 

290 

300 
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SUBHOUTIME MATPV ( A.B.C.NIU.NCA.NCB.irrA.mrB) 


trr 

IS 


m 

2S 

S(!NC 

** 

MT 

3S 

m MATRIX NULTIPLIGATION - GENERAL OR TRIANGULAR STORAGE 

** 

NTT 

40 



MT 

SO 

** CALLS SUBROUTINE LOC(SSP) 


MT 

66 

** 


MT 

70 

)K:k)K*)K)ltllt*)H*****»*)t^******)t‘************>li*H‘*************************** 

MT 

00 



MT 

90 

m'>» MATRIX IN GENERAL S10RAGB 


Mr 

100 

KT»1 MATRIX IN TRIANGULAR STORAGE (SYMMETRIC) 


MT 

110 

IMPLICIT BEAL»0(A-H,0-Z) 


Ml' 

120 

DIMENSION A( 1) . B(l) . C( 1) 


Ml' 

130 

DO 3 I> l.NRA 


MT 

140 

DO 2 K« 1 , NCB 


MI' 

130 

CALL LOG ( I.K, IT.HRA.NCB.O) 


MT 

160 

C( IT)»0.D0 


MT 

170 

DO 1 J»1,NCA 


Mr 

100 

CALL LOC ( I , J , I R, NBA, NCA, MTA) 


Ml' 

190 

CALI. LOC (J.K. IS.NCA.NCB.MTB) 


MI' 

200 

C< IT)«C( IT)+A< IR)*D( IS) 


MT 

210 

CONTINUE 


MT 

220 

CONTINUE 

CONTINUE 


MT 

230 


MT 

240 

RETURN 


MI' 

230 

END 


Ml' 

260 

SUBROUTINE LOC ( I . J , III, N, M.MS) 


LC 

10 


LC 

20 



LC 

30 

** sap SUBROUTINE - MATRIX STORAGE MANIPULATOR 

** 

LC 

40 

*!« 

** 

LC 

50 


LC 

60 



LC 

70 

IMPLICIT REAL*8(A-U,0-Z) 


LC 

00 

1X= 1 


LC 

90 

JX=J 


LC 

100 

IF(MS-l) l,2,S 


LC 

110 

TflX=N*(.IX“l) + IX 


LC 

120 

GO TO 7 


LC 

130 

IF( IX-JX) 3,4,4 


LC 

140 

IRX» IX+( JXifcJX- JX) /2 


LC 

130 

(X) TO 7 


LC 

160 

IHX=JX+( IX:«IX-IX)/2 


LC 

170 

GO TO 7 


LC 

100 

1RX=0 


LC 

190 

IF(IX-JX) 7,6,7 


LC 

200 

IRX=IX 


LC 

210 

IRa IRX 


LC 

220 

lUlTURN 


LC 

230 

END 


LC 

240 
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SOBROUTIIfB FRAME ( IX, I8X, lY, I8Y) 

FR 

10 


FR 

20 

FRAIKS A SCREEN WINDOW 

FR 

80 


FR 

40 

CALL NOVABS ( IX, IY> 

FR 

S0 

CALL DRWABS (IX, lY'flSYI 

FR 

60 

CALL DRWABS ( IX+ISX, lY'f ISY) 

FR 

70 

CALL DRWABS (IX+ISX, lY) 

FR 

00 

GALL DRWABS ( IX, IY> 

FR 

90 


FR 

100 

RETURN 

FR 

110 

END 

FR 

120 

SUBROUTINE UNITS (A,B,DX,DY) 

UN 

10 

UN 

20 

CONVERTS CENTIMETERS TO VIRTUAL COORDINATES 

UN 

30 


UN 

40 

IX«KCM<A> 

UN 

00 

IY»KCM(B) 

UN 

60 

DX«FLOAT< IX) 

UN 

70 

DY» FLOAT! lY) 

UN 

00 


UN 

90 

RETURN 

UN 

100 

END 

UN 

110 

SUBROUTINE RECT (DX,DY,D) 

Rd 

10 

Rd 

20 

DRAW A SOUARE WITH CENTER AT DX,DY 

Rd 

30 


Rd 40 
Rd 30 
Rd 60 
Rd 70 
Rd 80 
Rd 90 
Rd 100 
Rd 110 
Rd 120 
Rd 180 
Rd 140 
Rd 150 
Rd 160 
Rd 170 
EO 10 
Ett 20 


DRAW AN EQUILA'IERAL TRIANGLE 

WITH CENTROID AT DX.DY 

EO 

30 

CD - LENGTH OF TRIANGLE LEG 

EO 

40 



EO 

80 

PI»3. 14 


EO 

60 

DS»KCM( D> 


EO 

70 

CALL UNITS (DX,DY,S,T) 


EO 

00 

CALL MOVEA (S,T) 


EO 

90 

SSb . 433«DS 


EO 

100 

CALL MOVER (0. ,-SS) 


EO 

110 

CALL MOVER (-DS/2.,0.) 


EO 

120 

CAI.L DIUWR (DS.O.) 


EO 

130 

AW(ILE«PI/3. 


EO 

140 

K^COS( ANGLE) *DS 


EO 

180 

Y=S1W( ANGLE) :I!DS 


EO 

160 

G.\LL DRAWR (-X.Y) 


EO 

170 

CALL DRAim (-X.-Y) 


EO 

100 



EO 

190 

RETURN 


EO 

200 

END 


EO 

210 

SUBROUTINE CIRCLE (DX,DY.RS) 


CR 

10 



dll 

20 

DRAW A CIRCLE WI'IH CENTER AT 

DX,DY 

dll 

30 

R8 - RADIUS OF CIRCLE 

(Nl 

40 



dll 

80 


dD - LENCTO OF SIDE 
DS-KOHID) 

OALL WITS <DX,DY,8,T) 
dALL MOVEA (S.T) 

S«DS/2. 

CALL MOVER (-S.-S) 

CALL DRi\WH (DS.O.) 

CALL DRAWR (0.,DS) 

CALL DRAWR <-DS,0.) 

CALL DRAWR (0. ,-DS) 

RETCTWf 

END 

SUBROUTINE EOITTR (DX.DY.D) 
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Pl»8. 14 

CR 

60 

R»KCIf(RS) 

CR 

70 

CAUL UNITO (DX,DY,8,T) 

CR 

ao 

CALL flOVEA <S,T) 

CR 

90 

CALL NOVER (6..R) 

CR 

100 

C»2.*.PI*R+i 

CR 

no 

J«C 

CR 

lao 

AA»1./R 

CR 

100 

DO 1 l>l,J 

CR 

140 

A» I*AA 

CR 

ISO 

XxRiXBimA) 

CR 

160 

V»R!«COS(A) 

CR 

170 

CALL DBAWA (8+X,T^Y) 

CR 

100 

CONTINUE 

CR 

190 


CR 

200 

RETURN 

CR 

210 

END 

CR 

220 
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APPENDIX B 


VIP SAMPLE RUN 


B.l Introduction 

A VIP sample run Is presented In section B.2 to familiarize the 
user with the Interactive mode of operation and to Illustrate some of 
VIP's capabilities. This two-experiment session Is listed screen by 
screen (some are combined on one page to conserve paper) as viewed by 
the user and reproduced by the Tektronix hard copy unit. The entire 
session Is presented from LOGON to LOGOFF. Section B.3 Includes addi- 
tional output obtained by the user from the line printer (or VERSATEC) 
at the end of the session to serve as a record of that particular run. 

Experiments 1 and 2 address the question of observation corre- 
lations as discussed In section 3.2.6 and 4.3. As mentioned there. In 
an N-statlon configuration there are (N)(N-l)/2 baselines (and thus the 
same number of possible delay observations) but of those only N-1 Inde- 
pendent ones. In Experiment 1, a covariance analysis Is performed on a 
4-statlon network for the parameters described In section 3.2.3. Obser- 
vations from 3 Independent baselines are considered, their correlations 
determined according to the model of section 3.2.6. In Experiment 2, 
the same observation schedule Is followed but all 6 possible baselines 
are included. In this case, though, a diagonal observation covariance 
matrix is introduced (recall that applying the correlations between all 
six simultaneous observations would result In a singular observation 
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covariance matrix alnce the obaervatlons are linearly dependent). A 
comparison of the corresponding standard deviations of both experiments 
Indicates a decrease In those of Experiment 2 ranging from 20-30%. The 
largest differences Involve the earth orientation parameters since all 
baselines contribute to their estimation, In the global sense. It Is 
apparent that as the number of stations Increases and thus, the dis- 
crepancy between the total number of baselines and the number of inde- 
pendent baselines, the standard deviations become more optimistic when 
the correlations between simultaneous observations are neglected. 

To illustrate the effects of errors In the initial orientation 
of the baseline on the baseline components as explalm/d in Section 
3.2.3, perfect observations are simulated in each experiment. As can be 
seen in B.2, errors of 10 cm are introduced into the initial orientation 
of the pole and 1 ms of time in earth rotation over the first step. 

Those errors cause subsequent corrections to the approximate baseline 
’’components" of up to 29 cm in accordance with eqs. (3.2-24). 

The observation schedule and the simulated observations are 
given In B.3. The analysis considers a combination of delay and delay 
rate observations. The option of storing the observation schedule on 
file 9 (see Table A.l) prior to the session was chosen. The schedule 
was guided by two considerations: 

(1) that a source be observable simultaneously (maximum zenith distance 
of 80®) from all stations at a chosen epoch of observation, 

(2) that the final source schedule, over the 24-hour period of the 
simulations, be evenly distributed in right ascensions and declina- 
tions to achieve a strong geometry, especially to provide good 
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recovery for low and high aource declination-dependent parameters 
[Bock, 1980], 

In addition, the schedule Includes consideration of antenna slew time, 
cable wrap and tape constraints and was developed using a scheduling 
program written by Nancy Vandenberg at GSFC. The sources were selected 
from the source list obtained from the GSFC VLB! group and can be found 
in the sample run of B.2. 

At the end of Experiment 2, a typical visibility matrix is dis- 
played for the first five sources as viewed from the four participating 
stations. 
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*•2 A Typical Inf r«ctlv 8»»lon 

Th« following 1« a acreen-by-acraan display of Kicparlmanta 1 
and 2 from LOGOI? to LOGOFF. 
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iOXiPIATE STATION COORDINATES 
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B.3 Post-Session Output 


The following information is output (on file 7) by the line 
printer (or Versatec) as a record of a particular session. It can be 
collected by the user after the termination of the session. It con- 
tains all the output given in B.2 and in addition includes for each 
experiment : 

1. observation schedule (and optionally simulated observations) 

2. the normal matrix 

3. the variance-covariance matrix (unsealed) 

4 . the parameter correlation matrix 

5. the normal matrix eigenvalues in descending order. 

The observation schedule and simulated observations of Experi- 
ment 2 are presented on the following pages. The corresponding informa- 
tion for Experiment 1 is identical with those of baselines 1, 4 and 6 
of Experiment 2 at each epoch of observation. 

The parameter correlation matrix and normal matrix eigenvalues 
of Experiment 1 are also Included. Notice that the correlations are 
generally small indicating good separability of the parameters. The 
larg ?t correlations (0.8 -0.9) occur between the low-declination 
sources (also between each other) and the x(AX) component of the pri- 
marily east-west baseline. This follows from an examination of (3.3-4) 
(in this case the AY component of the Westford-Owens Valley baseline is 
relatively small) . 

The numbering of the parameters correspond to those of the standard 
deviations list given in B.2 
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The ratio of maximum/ iiiinimuin normal matrix eigenvalues is 
approximately 10® which is typical of multi-station configurations 
which are well-conditioned. 

Items 2 and 3 above are not presented here. 
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